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1.   INTRODUCTION  AND  PRELIMINARIES 

1.1  Introduction 

In  applications  such  as  reservoir  problems,  reactor  studies, 
and  numerical  weather  forecasting,  one  often  encounters  partial 
differential  equation  boundary  value  problems.  A  standard  technique 
for  solving  these  problems  is  to  approximate  the  solution  of  the 
differential  equation  at  a  discrete  set  of  points  by  the  solution  of 
a  linear  system,  Ax  =  b.   In  general,  such  systems  are  large,  sparse, 
and  often  nonsymmetric.   In  this  thesis  I  will  discuss  an  iterative 
method  for  solving  large  sparse  nonsymmetric  linear  systems. 

In  applications,  systems  arising  from  partial  differential 
equations  may  be  quite  large,  up  to  1,000,000  unknowns,  with  as  few  as 
5  to  10  nonzero  entries  per  equation.   The  size  and  spar si ty  of  these 
systems  motivate  the  use  of  an  iterative  method.   Direct  methods  require 
a  prohibitive  amount  of  storage. 

The  iteration  to  be  discussed  is  a  polynomial  based  gradient 
method  (Rutishauser  [7])  based  on  the  scaled  and  translated  Tchebychef 
polynomials .   In  the  following  I  will  show  that  when  applied  to  systems 
whose  spectrum  lies  in  the  right  half  complex  plane,  the  Tchebychef 
iteration  is  optimal,  in  a  certain  sense,  over  all  polynomial  based 
gradient  methods.   Further,  I  will  show  how  the  optimal  iteration 
parameters  associated  with  the  Tchebychef  iteration  can  be  computed 
from  knowledge  of  the  spectrum  of  the  system. 


A  drawback  to  most  iterative  methods  is  that  the  choice  of 
optimal  parameters  depends  upon  prior  knowledge  of  the  eigenvalue 
structure  of  the  system.   In  this  thesis  I  will  develop  an  adaptive 
procedure  for  determining  the  spectrum  of  the  matrix  during  execution 
at  a  relatively  low  cost.   Optimal  iteration  parameters  can  he  computed 
dynamically  with  little  prior  knowledge  of  the  spectrum. 

The  adpative  Tchebychef  algorithm  has  the  following 
properties: 

1.  Little  a  priori  knowledge  of  the  spectrum  is  required. 

2.  The  method  does  not  depend  upon  any  special  structure 
of  the  matrix  A. 

3.  The  method  is  sensitive  to  the  condition  of  the 

T 
matrix  A,  rather  than  the  matrix  A  A. 

h.      The  method  requires  only  one  matrix  vector  multiplication 

per  iterative  step. 

5.  The  method  can  be  used  in  conjunction  with  factorization 
methods . 

6.  The  method  may  be  used  on  singular  systems. 

Few  iterative  methods  have  been  developed  to  treat  nonsymmetric 
systems.  None  share  these  properties  with  the  adaptive  Tchebychef 
algorithm.   In  a  variety  of  test  problems,  the- Tchebychef  method  converged 
considerably  faster  than  two  competing  methods. 

Some  work  on  nonsymmetric  systems  has  been  done  by  Wrigley  [25], 
Kjellberg  [l6],  Kincaid  [15],  and  Faddeev  [8].   This  work  is  an  extension 
of  work  done  by  Diamond  [3]  and  Wrigley  [25].   Many  of  the  ideas  for 
this  thesis  came  from  work  done  by  Hestenes  [12],  Rutishauser  [7], 
Stiefel  [7],  [12],  [19],  Golub  [10],  and  Varga  [10],  [21]. 


The  main  ideas  of  this  thesis  are  as  follows:   Chapter  I 
presents  the  background  of  polynomial  based  gradient  metho 
Section  1.2  describes  the  type  of  system  to  be  considered.  A  crude 
region  containing  the  spectrum  of  the  system  is  determined.   Section  1.3 
outlines  the  general  polynomial  based  gradient  method  and  establishes 
criteria  for  choosing  a  sequence  of  polynomials  upon  which  to  base  an 
iterative  method. 

Chapter  2  deals  with  the  scaled  and  translated  Tchebychef 
polynomials  in  the  complex  plane.   Section  2.1  describes  the  properties  of 
the  scaled  and  translated  Tchebychef  polynomials.   It  is  shown  that  the 
asymptotic  behavior  of  the  Tchebychef  polynomials  is  related  to  ellipses 
in  the  complex  plane.   Sections  2.2,  2. 3,  and  2.h   show  that  the  scaled 
and  translated  Tchebychef  polynomials  fit  the  criteria  established  in 
Section  1.3 «   The  iteration  based  on  the  Tchebychef  polynomials  is 
described  in  Section  2.k   where  it  is  shown  that  the  iteration  parameters 
can  be  determined  from  the  scaling  and  translating  parameters  c  and  d. 

The  problem  of  choosing  the  parameters  c  and  d  is  outlined  in 
Chapter  J.   In  Section  3.1  a  measure  of  the  rate  of  convergence  is 
established  at  each  eigenvalue  as  a  function  of  the  parameters  c  and  d. 
A  mini -max  problem  is  constructed  in  terms  of  the  eigenvalues  and  the 
parameters  c  and  d  whose  solution  represents  the  "best"  choice  of 
parameters  c  and  d.   The  mini -max  problem  is  refined  and  reformulated 
in  Section  3.2.   In  Section  3-3  it  is  shown  that  for  a  real  matrix  A  the 
iteration  can  be  carried  out  in  real  arithmetic. 


Chapter  k   is  devoted  to  finding  the  solution  to  the  mini -max 
problem  constructed  in  Chapter  3'   Section  k.l   shows  that  the  solution 
of  the  mini-max  problem  can  be  found  in  terms  of  the  solution  of  a 
mini -max  problem  when  the  spectrum  of  the  matrix  A  contains  1,  2,    or 
3  complex  conjugate  pairs  of  eigenvalues.   The  remainder  of  the  chapter 
is  devoted  to  solving  the  mini-max  problem  in  these  special  cases. 
Section  k.k   contains  the  algorithm  for  solving  the  mini-max  problem. 

Chapter  5  describes  a  method  of  estimating  the  eigenvalue 
structure  of  the  matrix  A  during  execution  and  a  dynamic  procedure  for 
improving  the  choice  of  parameters  c  and  d.   Section  5-1  deals 
with  an  adaptation  of  the  modified  power  method  for  simultaneous 
estimation  of  several  eigenvalues  based  on  the  sequence  of  residuals 
(Wilkinson  [2^]).   Section  5-2  shows  how  each  new  eigenvalue  estimate 
can  be  added  to  previous  estimates  to  yield  an  improved  choice  of 
parameters  c  and  d. 

Chapter  6  gives  a  discussion  of  the  implementation  of  the 
algorithm  followed  by  a  discussion  of  competing  methods  and  experimental 
results.   It  is  shown  that  the  adaptive  Tchebychef  algorithm  requires 
fewer  iterative  steps  to  produce  convergence  and  half  as  much  work  per 
step  than  two  competing  methods:  Bidiagonalization  (Golub  and  Kahan  [11]) 

and  the  method  of  Conjugate  Gradients  (Hestenes  and  Stiefel  [12]) 

T      T 
applied  to  the  equivalent  system  A  Ax  =  A  b. 

The  appendix  contains  a  listing  of  the  algorithm,  coded  in 
FORTRAN . 


1.2  Type  of  System 

The  adaptive  Tchebychef  algorithm  is  restricted  to  system 
with  eigenvalues  in  the  right  half  (left  half)  complex  plane.   Such 
systems  arise  naturally  from  applications  to  partial  differential 
equations,  especially  in  conjunction  with  finite  element  methods 
(Zienkiewicz  [27])  . 

In  the  remainder  of  this  thesis  the  matrix  A  will  be  a  real 
valued  matrix  with  eigenvalues,  X. ,  in  the  open  right  half  plane. 
The  possibility  of  zero  eigenvalues  will  be  discussed  in  Chapter  6. 

To  establish  bounds  on  the  spectrum  of  such  a  matrix,  let 

T  T 

A+A         A-A 
M  =  — 5 —  and  N  =  — - —  .   It  is  clear  that  M  is  symmetric,  N  is  anti- 
symmetric, and  A  =  M  +  N.   Suppose  \   is  an  eigenvalue  of  A  corresponding 
to  the  eigenvector  v  =  x  +  iy.  Then, 

<AV,V>   <MV,V>   <NV,  V> 
<v,v>  ~~  <v,  v>  +  <v,v>  ' 

Because  M  is  symmetric,  we  have 

<Mv,v>  =  <M(x+iy),  (x+iy)>  =  <Mx,  x>  +  i<My,  x>  -i<Mx,y>  +  <My,y> 
=  <Mx,x>  +  <My,y>. 

Since  N  is  antisymmetric,  <Nx, x>  =  <Ny,y>  =  0,  and  we  have 

<NV,  v>  =  <N(x+iy),  (x+iy)>  =  <Nx,x>  +  i<Ny,  x>  -  i<Nx, y>  +  <Ny,y> 
=  2i<Ny,x>. 

Thus,  we  have 


Re(?0  3  <**,*>  +  <My,y>  } 


and 


|lm(\) 


2  [<Ny,  x>  | 
|x||2  +  ||y||2 


< 


2   N 


x 


|y| 


|x||2  +  ||y||2      " 


<  N 


Since  M  is  symmetric,  its  spectrum  is  contained  in  some 
interval,  [a,b],  on  the  real  line.   Since  N  is  antisymmetric,  its 
spectrum  lies  along  the  imaginary  axis  and  ||n||  is  equal  to  the  spectral 
radius  of  N.   If  \  is  an  eigenvalue  of  A,  then 

Re(\)  e  [a  b]  , 
|lm(X)  |  <  ||n||  . 

Using  Gershgorin's  circle  theorem,  bounds  for  the  interval  [a,b]  and  ||n|| 
can  be  found  in  terms  of  the  elements  of  the  matrices  M  and  N  (Varga  [22]) 
giving  a  rectangle  known  to  contain  the  spectrum  of  the  operator  A  (see 
Figure  1.1  and  Figure  1.2). 


Figure  1.1 


Figure  1 . 2 


Not      at  if       small  with  respect  to  |b-a|,  the 
rectangle  is  like  that  shown  in  Figure  1.1.   If  ||n!!  is  large  with 
respect  to  |b-a|,  then  the  rectangle  is  like  that  shown  in  Figure  1.  ■ 
Both  rectangles  are  more  closely  approximated  by  an  ellipse  than  a 
circle.  The  importance  of  this  observation  will  be  shown  later. 

If  M  is  positive  definite,  then  the  spectrum  of  A  lies 
the  right  half  plane.  In  particular,  if  A  is  diagonally  dominant  with 
positive  real  diagonal  elements,  then  M  will  be  positive  .1  ;i. .  - 
This  type  of  matrix  is  often  encountered  in  application.  Of  r     _ar 
interest  are  positive  definite  systems  perturbed  by  a  small  nonsynmetric 
part,  yielding  a  system  with  its  spectrum  contained  in  a  region  like 
that  shown  in  Figure  1.1.   It  is  this  type  of  system  for  which  the 
Pchebychef  algorithm  has  the  greatest  advantage  over  competing  methods . 

1.;  Iype  *>f  Iteration 

hyp;;;   ;   ;.:.-  :e  solve  the  system  Ax  =  "b.   If  x_  is  a:: 
initial  guess  at  the  solution,  x,  we  car.  .-er.pvte  the  initial  res 

rQ  =  A(x-x^  =  c  -Ax.  . 

and  use  this  information  to  make  a  new  guess,  >:.  .  Again  we  ear.  :smpute 
the  residual 

r  =  A(x-x  )  =  b  -Ax 

1-1 

Utilizing  all  of  the  information  at  hand  we  can  use  both  r  and  r  to 

3  a  new  guess,  x„ .   Ihis  leads  to  an  iteration  with  general  step 


8 


n 

x  ,  =  x  +  Z   7  .  r.  . 
n+1    n   .    'm  1 


where  the  y.    ' s  are  constants.  At  each  step  we  add  a  linear  combination 
of  the  previous  residuals.   Let  e  =  x-xn  be  the  error  at  the  nth  step. 


From  the  general  step  we  can  write 


n  n 

e  -,  =  e  -  Z  7  .  r.  =  e  -  A  2  y   .  e.  . 

n+1    n   .  ,  'm  i    n     .  _  'm  i 
i=l  i=l 

Lemma  1.1  The  error  at  the  nth  step  is  given  by 

e  =  P  (A)en  , 
n    n    0  ' 

where  P  (z)  is  a  polynomial  of  degree  n  such  that  P „(0)  =  1. 
Proof  The  proof  will  proceed  by  induction.   For  n  =  1  we  have 

ei=  eo  "  ?oo  Aeo  =  (I"7oo  A)eo  ' 

so  that  we  may  define 

Px(z)  =  (l-700  z)    , 

and  we  have  P  (0)  =  1. 

Assume  the  conclusion  is  true  for  i  <  n;  then, 

n 
e   -,=e   -A  y      y    .    e. 
n+1    n     *  'm  i 
i=l 

n-1 
=  e  -  y       Ae-AZy.e. 
n   'nn   n     ._  'm  1 

=  [(1-7   A)P  (A)  -  A  ^  7  •  ?.  (A)l  e   . 
nn    n        .  .,   ni  -1 
i=l 


n-1 

P  _(z)  =  (1  -7   z)P  (z)  -  z  L     7    .    P.(z)   . 
n+lv         nn    n       .  ,   m  i 

1=1 


Then,  P   (z)  is  a  polynomial  of  degree  n+1  and  P   (0)  =  1.  This 
proves  the  lemma. 

Any  sequence  of  polynomials,  {P  (z)},  can  be  generated  by 

choosing  the  constants,  [7. .}.  We  would  like  to  choose  the  sequence 

J-  J 


of  polynomials  so  that 


Henl|<||Pn(A)||    ||e0H 

is  small.  Let  us  examine  P  (A).  If  A  is  diagonalizable,  then  the 
Jordan  form  of  A  is  a  diagonal  matrix  (Birkhoff  and  MacLane  [l]). 
There  exists  a  nonsingular  S  such  that  A  =  S~  JS  and 


\ 


Now 


J  = 


P  (A)  =  P  (S_1JS)  =  S'1?    (J)S  , 


and  since  j  is  a  diagonal  matrix  we  have 


Pn(j) 


This  yields  the  following  result: 

Theorem  1.2     If  A  is  diagonalizable,    then   ||P   (A)  ||  -»  0  as  n  -♦  co  if  and  only 

if  P  (\.)   -*  0  as  n  -*  00   for  every  eigenvalue,    \. ,    of  A. 
n     1  '      1* 


10 


Proof  The  proof  is  clear  from  the  discussion  above. 


Suppose  A  is  not  diagonalizable;  that  is,  suppose  A  has  non- 
linear elementary  divisors.   The  Jordan  form  of  A  has  nontrivial 
Jordan  blocks .  We  have  A  =  S~  JS  where 


Ji 


j 


Jk 


and 


X.        1 

1 


J.  = 
l 


is  the  Jordan  block  associated  with  the  eigenvalue  X.    of  multiplicity 

d. .   In  this  case 

i 


P  (A)  =  P  (S_1JS)  =  S_1P  (J)S  , 
n      n  n     ' 


and 


P  (J)  = 
n 


\ 


P  (JJ 

n  k 


Taking  successive  powers  of  J.  we  see  that 


2 
\.  Zk.  1 

i       l 


j2- 


J?- 

1 


IK 


and  in  general, 


cX"1    (>r2 


(s   )xs-di+l 
di     x 


v2;    l 


•  (Xs-1 


l 


s  .  8-1       s(s-l)  ,  s-2  s(s-l)  . .  .  (s-r+1)   .  s-d-;+l 


s(s-l)%s-2 
2T~\ 


s  .  s-1 


12 


The    superdiagonal  terms  act  like  derivatives.      Since  P   ( J. ) 

n     1 

is  a  linear  combination  of  powers  of  J.,    we  have 


'p  (x.)     p«(x.)     5tP"(\.)    .  .  .     /.  1lU  p ri'^oo 

n\   x;  nv   i;        2!    nv   i;  (d.-l)!    n  v    i 


P   (J.)= 

nv   i7 


This  yields  the  following  theorem: 
Theorem  1 .3  If  V-  is  a^  eigenvalue  of  A  with  multiplicity  d.,  then 
|jPn(A)  ||  ->  0  as  n  -  oo  if  an<i  only  if  P^  (\. )   -*   0  as  n  -*  oo  for  every  j  <  d., 
for  each  eigenvalue  A..  . 
Proof  The  proof  is  clear  from  the  discussion  above. 

When  looking  for  a  sequence  of  polynomials  to  suppress  the 
eigenvalues  of  A,  we  must  find  one  whose  derivatives  also  suppress  the 
eigenvalues  of  A. 

In  light  of  the  previous  discussion,  we  can  establish  three 
criteria  upon  which  to  choose  a  sequence  of  polynomials: 

1.  We  must  choose  P  (z) ,    among  polynomials  of  like  degree 
such  that  P  (0)  =  1,    to  be  "as  small  as  possible" 

on  the  spectrum  of  A. 

2.  If  A  has  nonlinear  elementary  divisers,  then  we  must 
choose  Pn(z)  so  that  its  derivatives  are  small  on  the 
spectrum  of  A. 


}•  We  must  choose  P  (z)  to  have  some  recursive 

properties  so  that  all  of  the  previous  residuals 
need  not  be  stored. 

In  the  next  chapter  we  will  see  that  the  scaled  and  translated 
Tchebychef  polynomials  fit  these  criteria. 


Ik 


2.   TCHEBYCHEF  ITERATION  IN  THE  COMPLEX  PLANE 

The  Tchebychef  polynomials  were  discovered  a  century  ago  by 
the  Russian  mathematician  Tchebychef  (the  spelling  of  which  has  many 
variations) .   Their  importance  for  practical  computation,  however,  was 
rediscovered  about  thirty  years  ago  by  C.  Lanczos.   Since  then  they  have 
found  many  uses  in  numerical  analysis  (Fox  [9l)« 

The  definition  and  basic  properties  of  the  Tchebychef  polynomials 
in  the  complex  plane  will  be  discussed  in  Section  2.1.   Sections  2.2, 
2.3,  and  2.k   will  show  how  the  Tchebychef  polynomials  meet  the  criteria 
established  in  Section  1.3«   The  gradient  method  based  on  the  Tchebychef 
polynomials  is  developed  in  Section  2.k. 

2.1  The  Tchebychef  Polynomials 

The  Tchebychef  polynomials  are  given  by: 

T0(z)  =  1, 


Tx(z)  =   z, 
Tn+1(z)  =  2zTn(z)  -  Tn-1(z),  n  >  1 


They  may  also  be  written: 

T  (z)  =  cosh(n  cosh"  (z)), 

where  the  branch  of  cosh"  with  positive  real  part  is  used. 


L5 


Consider  the  map  t]  =  cosh(z)  .      Let   z  =  x  t  iy,    T]        u   t  Lv. 
Then  cosh(z)   =   cosh(x+iy)   =  u  +  iv  =  r\ .      U  "ing   the  expansion    Cormula 
for  the  cosh,    we  have 

cosh(x+iy)   =   cosh(x)    cos(y)    +   i   sinh(x)    sin(y), 


or 


u  =   cosh(x)   cos(y), 
v  =   sinh(x)    sin(y) . 

Suppose  we  fix  x  >  0  and  allow  y  to  vary.   Then  u  and  v 


satisfy 


2         2 
u         v 

+  -z =  1  • 


2         2 
cosh  (x)    sinh  (x) 


That  is,  the  line  x  =  constant  is  mapped  onto  an  ellipse  with  semi- 
major  axis  |cosh(x)|,  semi-minor  axis  |sinh(x)|,  and  foci  at  1  and  -1 
(see  Figure  2.1).   This  map  has  period  2rti.   Since  cosh(x)  and  sinh(x) 
are  increasing  for  x  >  0,  if  0  <  a  <  b,  then  the  line  x  =  a  is  mapped 
onto  an  ellipse  inside  and  confocal  to  the  ellipse  that  the  line  x  =  b 
is  mapped  onto  (see  Figure  2.1).   If  x  =  0  we  have 

u  =  cos(y), 

v  =  0  , 

and  the  ellipse  has  collapsed  onto  the  real  line  segment  [-1,  1] . 

Because  of  periodicity,  the  map  r\   =  cosh(z)  takes  the  region 

shown  in  Figure  2.2  onto  the  entire  t}- plane.   Each  vertical  line  in  this 

region  is  mapped  onto  an  ellipse  in  the  t} -plane.   This  region  is  the  branch 

of  cosh"  used  in  the  definition  of  the  Tchebychef  polynomials. 

The  function  cosh"^  may  also  be  written  in  log  form: 

1 
cosh"  (w)  =  ^ (w  + (w  -l)  )  . 
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COSH(b) 


x=a        x=b 


Z  -  PLANE 


COShfV/j) 


COSH(z) 


rj  -PLANE 


Figure  2 . 1 


2  7TI  -- 


7TI 


Figure  2.2 


re  must  be  taken  when  choosing  the  branch  of  the  square  root.  The 
branch  chosen  depends  on  the  argument  w  and  should  be  cho  i  that 
(w  )tl  =   w. 

The  nth  Tchebychef  polynomial,  T  (z)  =  cosh(n  cosh"  (z)), 
maps  an  ellipse  onto  a  vertical  line  segment  in  the  region  shown  in 
Figure  2.2,  multiplies  this  line  segment  by  n,  and  maps  the  new  line 
segment  back  onto  another  ellipse  (see  Figure  2.1).   Since  the  new  line 
segment  cuts  through  n  branches  of  cosh  ,  it  is  wrapped  around  the  new 
ellipse  n  times.  The  degenerate  ellipse,  the  line  segment  [-1,  1],  is 
mapped  onto  the  line  segment  [0,  in],  multiplied  by  n,  and  mapped  back 
onto  the  line  segment  [-1,  l] .   Since  it  is  wrapped  around  n  times, 
Tn(z)  has  n  zeros  on  the  line  segment  [-1,  l] . 

To  establish  some  notation,  let  ^(d,  c)  be  the  family  of 
ellipses  centered  at  d  with  foci  at  d+c  and  d-c.   Let  F(d,  c)  e  <?(d,  c) 
be  a  member  of  this  family.   Let  F.(d, c)  c  F . (d, c)  mean  that  the  ellipse 

X  J 

F.(d, c)  is  inside  the  ellipse  F.(d, c).   Let  z  e  F. (d,  c)  mean  the  point  z 
is  on  the  ellipse  F. (d,  c).   The  Tchebychef  polynomials  then  map  members 
of  ^(0,1)  onto  other  members  of  7^0,1). 
Lemma  2.1  Suppose  z.  eF.(0,  l),  z.  eF.(0,  l);  then, 

A     J.  J     J 

Re(cosh_1(z.))  <  Re  (cosh'1  (z.))  OF.  (0,1)  c  F.(0,l), 

Re(cosh_1(z.))  =  Re(cosh_1(z.))  ^>  F.  (0,1)  =  F.(0,l). 
■*■  J  J 

Proof  The  proof  follows  from  the  discussion  above. 

Now  consider  the  scaled  Tchebychef  polynomials, 


T  (z) 


nv  0 
These  polynomials  exhibit  an  asymptotic  behavior 
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Tn(z) 
Lemma  2.2     If  C    (z)    =  -^—, r  ,     z_   I  [-1,    l],    then  for  large  n 


nv   0' 

C    (z)   =  en^cosh~    ^  -cosh"    (zo)l 
n^ 


Proof     From  the  definition  of  the  cosh,  we  have 


n  cosh"    (z)  -n  cosh"    (z) 

Cn(z)   =  " 1 — 1 

gn  cosh     (zQ)+e-n  cosh     (zQ) 


Since  the  branch  of  cosh-1  with  positive  real  part  is  used  the  result 

follows . 

For  large  n,    C    (z)    takes  on  the   form  r  ,    which  motivates  the 

following  definition. 

Definition     Let 

1 

r(z)   =  lim    |Cn(z)n| 

Z-+  °° 

be  the  asymptotic  convergence  factor  of  C  (z)  at  the  point  z. 

The  asymptotic  convergence  factor,  which  will  be  referred  to 
as  the  convergence  factor,  is  related  to  the  rate  of  convergence  as 
defined  by  Young  [26]  and  Varga  [22]  in  that  rate  of  convergence  equals 

-  &n(r(z))  .      From  the  lemma  above  we  have 

r(z)  =  |eCOsh~  ^  "cosh"  (zo">  I  =  eRe(cosh"  (z))  -Re(cosh"  (zQ)) 

There  is  a  relationship  between  the  convergence  factor  and  the  members 
of  #(0,1),  the  family  of  ellipses  with  foci  at  1  and  -1. 

Lemma  2.3  If  z_  eF_(0,l),  z.  eF.  (0,1),  z.eF.(0,l),  then 

0   0V  '    "      l   i  '  ■  •  J    3 

r(z.)  <  r(Zj.)  &  F.(0,1)  c  F^(0,1), 


r(z1)  =  r(z.)  <>  P1(i  ,!)   K.(0,1), 
r(z)  =  1  O  z  e  FQ(0,1)  . 

t^   -  o-     t   \  Re  (cosh"  {?,))   -Re  (cosh"  (zn))   ..      ..  _  . . 

Proof  Since  r(z)=e         v  "  ^-U",  the  result  follows 

from  the  previous  lemmas. 

Tn(z) 
Theorem  2.U  Let  CQ(z)  =  T  ^  y,  zQ  /  [-1,  1] .   If  FQ(0,1)  is  the 

member  of  the  family  9^(0,1)  passing  through  z  ,    then 

f   0    if  z  is  inside  F_(0,l) 
lim  Cn(z)  =  ° 

n-»oo  l^  co    if  z  is  outside  F  (0, 1)  . 

Proof  From  Lemma  2.3  we  see  that  r(z)  <  1  if  z  is  inside  F  (0, 1)  and 
r(z)  >  1  if  z  is  outside  F  (0,1).   Since  |c  (z) |  =  r(z)n  the  result 
follows. 

Consider  the  transformation  z  =  -=—,    z^  =  — ,  where  d  and  c 

c  '      0   c 

are  any  complex  numbers.   Let 

T  &h) 

P  (X)  =  C  (z)  =  n  g    . 

n       n       T  (-) 
n  c 


The  polynomials  P  (\)  are  called  the  scaled  and  translated  Tchebychef 
polynomials.  Notice  that  P  (0)  =  1  as  is  required  for  use  in  a  gradient 
method. 

The  choice  of  d  and  c  determines  a  family  of  ellipses,  <£*(d,  c), 
with  foci  at  d+c  and  d-c  (see  Figure  2.3).   Let  F„(d,c)  e  *^(d,c)  be  the 
member  of  the  family  passing  through  the  origin.   The  transformation 
maps  members  of  ^(d,  c)  in  the  \-plane  onto  members  of  ^(0,1)  in  the 
z-plane . 
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F0(^c) 


Fo(0J) 


X-PLANE 


Z -PLANE 


Figure  2-3 


As  before,    let 


r(\)   = 


lim    |P  (\) 
n-*oo 


be  the  asymptotic   convergence   factor  of  P   (\)    at  the  point  \.     We  have 


from  above 


,.v  Re(cosh'1(— ))  -ReCcbsh"1^)) 

r(\)   =   e      v  v    c  c 


The  relationship  between  the  members  of  *&(&,  c)  and  the 


convergence  properties  of  the  polynomials  P  (\)   is  given  in  Theorem  2.5 

t  {£±) 

Theorem  2.^     If  P   (X)   =  — — 2 — ,    then 

n  T   (S) 

nvc 


.  ! 


1.  f  0   if  X  ia  inside  F  (d,  c) 
Li*  V  (\)   = 

n-  ^  *   if  X  is  outside  F  (d,c) 

If  X.  eF.(d,c),  X.  €F.(d,c),  then 
l   iv  '  '   j   j 

2.  r(\.)  <  r(X  )  <>  Fi(d,c)  c  Fj(d,c), 
r(\i)  =  r(Xj)  o  F.(d,c)  =  Pj(d,c), 
r(X)  =1     <=>  X  €  FQ(d,c)  . 

Proof  Since  the  transformation  maps  members  of  ^(d,  c)  in  the  X-plane 
onto  members  of  ^(0,1)  in  the  z-plane,  the  result  follows  from 
Theorem  2.k   and  Lemma  2.3 • 


2.2  Optimal  Properties  of  the  Tchebychef  Polynomials 

The  first  criterion  mentioned  in  Section  1.3  suggests  that  when 
choosing  a  sequence  of  polynomials  upon  which  to  base  an  iterative 
method,  it  is  desirable  to  choose  polynomials  that  are  small  on  the 
spectrum  of  the  matrix  A.   Since  the  spectrum  of  A  is  seldom  known, 
it  is  more  practical  to  choose  polynomials  that  are  small  on  a  region 
containing  the  spectrum  of  A.   If  the  region  is  bounded  by  a  circle  or 
an  ellipse,  the  scaled  and  translated  Tchebychef  polynomials  have  certain 
optimal  properties.  Much  is  known  of  the  optimal  properties  of  the  monic 
scaled  and  translated  Tchebychef  polynomials  over  all  monic  polynomials 
of  like  degree  on  regions  bounded  by  ellipses  (Hi lie  [13],  Walsh  [23]). 
Similar  results,  but  not  as  strong,  can  be  shown  for  polynomials 
normalized  at  the  origin. 

Definition  Let  S  =  {all  polynomials,  s  (X),  of  degree  n  such  that 

s  (0)  =  1} .   The  elements  of  S  are  said  to  be  normalized  at  the  origin. 

nv  '  '  n — 
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Theorem  2.6  Let  E  be  a  closed  and  "bounded  infinite  set  in  the 
complex  plane.   There  exists  a  unique  t  eS  such  that 

max  |t  (z) |  =  min   max   |s  (z) | . 

zeE  s  eS    Z€E 

n  n 

Proof  Omitted  (Hi lie  [13])- 

If  the  region  E  is  bounded  by  an  ellipse,  a  circle  being  a 
special  case  of  an  ellipse,  then  the  maximum  will  occur  on  the  boundary. 
Using  the  notation  of  Section  2.1,  let  F  (d,  c)  be  the  member  of  the 
family  ^(d,  c)  with  semi -major  axis  a  <  0.   Instead  of  taking  the  maximum 
over  the  entire  region  we  may  take  the  maximum  over  the  boundary, 
F  (d,  c).   The  circle  with  center  d  and  radius  a  is  denoted  F  (d, 0). 
If  the  spectrum  of  A  is  contained  in  a  region  that  is  bounded  by  a  circle 
that  does  not  include  the  origin  in  its  interior,  we  have  the  following 
result. 

Theorem  2.7  Suppose  F  (d, 0)  does  not  include  the  origin  in  its  interior. 
Let  a  <  |d|.   The  unique  polynomial  tn  e Sn  such  that 


max    |t  (X)  j  =  min      max    |s  (X) 
XeF   (d,0)   n       s  eS   XeF   (d,0)   n 


a     '  n     n  a 


is  given  by 


y*>  -  (^)n 


Proof     Suppose   q     eS     and 

max        |q    (\)  |   <         max        |t    (X)  |    =    (-4T)n 
7   (an)    ^  x?v  (c\.o\     n  \lal/ 


\€Fa(d,0)  ^Fa(d,0) 


Then, 

|q(\)|    <    |t(X)|    -      J^jn    , 

for  every  X  eF   (d,  0) .      By  Rouch^'s  theorem   the  polynomial  t    (X)  -  a(X) 
has  the   same  number  of  zeros  inside  F   (d,  0)   as  the  polynomial 

8 

t  (\)  does.  Notice  that  t  (0)  -a  (0)  =  0.   Since  X   =  0  is  not  in 
the  interior  of  F  (d, 0),  t  (X)   -q  (X)  is  a  polynomial  of  degree  n  with 
n+1  roots.  We  can  conclude  that  t  (\)  =  q  (X) ,    and  the  theorem  is 
proved. 

If  the  region  is  bounded  by  an  ellipse  with  real  foci  that 
does  not  contain  the  origin  in  its  interior,  we  have  the  following 
result  due  to  Clayton  (Wrigley  [25]). 

Theorem  2.8  Let  0  <  c  <  a  <  d.   The  unique  polynomial  t  eS  such  that 


max   |t  (X)  |  =  min     max   |s  (X) 

\eF  (d,c)   n       s  eS  XeF   (d.c)   n 

a  n  n     av  '    ' 


is  given  by 


T  (^) 
t  (X)   =  P  (\)  =  n  °   , 

n        n        T  (1) 
n^c 


the  associated  scaled  and  translated  Tchebychef  polynomial. 
Proof  emitted  (Wrigley  [25 ] ) . 

This  result  cannot  be  extended  to  d  and  c  with  complex  values 
(an  easy  example  can  be  constructed),  but  it  can  be  shown  to  be 
asymptotically  true.   For  large  n,  P  (X)   tends  very  quickly  to  the 
optimal  polynomial  in  S  . 


2k 


Lemma  2.9  Suppose  F  (d, c)  does  not  contain  the  origin  in  its  interior. 

■  3. 

Let  t  e  S  be  the  unique  polynomial  such  that 

max   |t  (X)  |  =  min     max   |s  (\) | ; 
XeFQ(d,c)   n       seS^  XeFo(d,c)   n 


n  n     a 


then, 


min   |P(>0|<    max   |t  (X)  |  <    max   |p(X)  | 
XeF  (d,c)  XeF  (d,c)   n       XeF  (d,c) 

a.  a  d 


Proof  The  second  inequality  is  true  by  hypothesis.   Suppose  that 


max   |t  (X)  I  <    min   |P  (X) |; 
XeFa(d,c)   n      XeFa(d,c) 


then, 


t  (X)  <  P  (X) 

nv      nv  ' 


for  every  X  eFa(d,c).   By  Rouche's  theorem,  the  polynomial  P  (X)  -t  (X) 
a  n     n 

has  the  same  number  of  zeros  in  the  interior  of  Fa(d,  c)  as  Pn(X)  does. 

P  (X)  has  n  zeros  on  the  line  segment  joining  the  foci,  d+c  and  d-c. 

Notice  that  P  (0)  -t  (0)  =  0.   Since  X  =  0  is  not  in  the  interior  of 

F  (d, c),  P  (X)  -t  (X)  is  a  polynomial  of  degree  n  with  n+1  zeros.  We 
a       n     n 

can  conclude  that  P  (X)  =  t  (X),  and  the  lemma  is  proved. 

Theorem  2.10  Suppose  F  (d, c)  does  not  include  the  origin  in  its  interior. 
Let  t  e  S  be  the  unique  polynomial  such  that 

max   (t  (X))  =  min     max   (s  (X))  . 

XeF  (d,c)   n       s  €S   XeF  (d,  c)  n 
av  '  n  n     av  ' 
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Let 


M(s  )  ■    max    lsn(^  I  ; 
^€F.(d,c) 


then, 


lim  [M(tn)n]  =  lim  [M(Pn)n] 
n-»  oo  n-+oo 


Proof  Let 

m(P  )  =    min   |p  (\) |. 
n    \€Fa(d,c)  n 

From  Lemma  2.9  we  know  that  m(P  )  <  M(t  )  <  M(P  ).   It  is  sufficient 

n  —    n  —    n 

to  show  that 

1  1 

lim  [m(P  )n]  =  lim  [M(P  )n]  . 
n  n 


By  Theorem  2.5  all  points  on  the  ellipse  F  (d, c)  have  the  same  convergence 

factor;  thus,  we  have 

1  1 

r(\)  =  lim  [m(Pn)n]  =  lim  [M(Pn)n] 
n.-+oo  n-*  °o 


for  every  X  eF   (d, c).   This  proves  the  theorem, 
a 

Because  of  the  nature  of  the  cosh  function,  the  asymptotic 
convergence  factor  is  achieved  very  quickly;  thus,  the  scaled  and 
translated  Tchebychef  polynomials  tend  very  quickly  to  the  optimal 
polynomial  in  S  . 

As  the  focal  length  c  approaches  0  the  ellipse  F  (d, c)  is 
deformed  into  the  circle  F  (d, 0).   To  show  that  the  result  for  circles  is 

3. 

compatible  with  the  result  for  ellipses  we  have  the  following  lemma: 
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Lemma  2.11  If  P  (\)  =  — — § —  ,  then 

n       T  £)■ 

n  c 


lim  P  (X)  =  (^)n 
c->0 


Proof  By  the  definition  of  cosh,  we  have 


T  (±±)  n  cosh"1^)   -n  cosh"1^) 

P  (x)  =    n    ?  =  * c    +  e L- 

n      T  (-)  -lA         -lA 

nvc7      n  cosh   ("cJ    -n  cosh  {■$) 


which  tends  to 


-1  — 
n  cosh  (  c  > 

n  cosh"  (77) 


-1      -1  2   2 

as  c  decreases.   By  the  log  form  of  cosh   (cosh"  (wl  =  &?(w  + (w  -1)  ))  this 


is 


c        c  

d/c  +  [(d/c)2  -  1]1/2 


(d-\)  +  [(a-\f  -  c2]l/2 

d  +  [d  -  c  J  ' 


-d-\ 


which  tends  to  (— jr-)   as  c  decreases.   This  proves  the  lemma, 


2-3  Convergence  of  P^(\) 
t>        n 

Recall  from  Section  1.3  that  if  the  matrix  A  has  nonlinear 
elementary  divisors,  the  derivatives  of  the  sequence  of  polynomials  must 
also  converge  to  zero  on  the  eigenvalues  of  A.   In  this  section  it  will 
be  shown  that  each  derivative  sequence  of  P  (\)  has  the  same  region  of 
convergence  as  the  sequence  P  (\)   does. 
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T  ( ) 

Theorem  2.12  Let  P  (\)  ■  — — ^ —  .   [f  \   is  inside  V   (d, c),  bhe 

D       T  (-) 

n  c 

member  of  ^(d,  c)  passing  through  the  origin,  then 


lim  P^(\)  =  0  , 


for  every  j  . 

Proof     Suppose  \  is  inside  F   (d,  c).      Since  the   interior  of  F„(d,  c)    is 
an  open  set,    there  is   some  5  >  0  such  that  T   =    {t/|t-\|   =   8}    is  inside 


Fn(d,  c)  .     We  have 


v^(\)   =  -^-     f        ?n  dt 


r 


Let  |P  (\  ) |  =  max   |p  (\) |;  then, 


n  n     \€T    n 


t  -\  ..    r    p  (t) 


<31  'vv  rdt 


|P  (X    ) 

.,  '  n  n 


(5)' 


Since  r  is  in  the  interior  of  Fn(d,  c)  we  know  that 


|P  (\  )  I  =  rn 
1  nv  ny  ' 


for  some  r  <  1,  and  the  theorem  is  proved. 
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The    speed  of  convergence  of    |p   (X  ) |    depends  upon  how  close 
T   is  to   the  boundary,    F   (d,  c).      There  is  a  trade  off  between  choosing 
smaller  5  and  the   speed  at  which    |p   (X   )  |    converges.      One  can,    in  fact, 
pick  a  new  5   for  each  n  and  get 

IP^OOI   <K(\,J)  nJ  |Pn(\)| 

where  K(A.,  j)  is  a  constant  depending  on  X   and  j.   From  Section  2.1  we 
know  that  |p  (A.)  |  =  r(X)   ,    so  that 

|p£°(X)|  <K(\,d)  nJr(\)n  . 

We  can  conclude  that  when  the  ellipse  F  (d, c)  contains  the 
spectrum  of  the  matrix  A  in  its  interior,  a  gradient  method  based  upon 
the  associated  scaled  and  translated  Tchebychef  polynomials  will 
converge,  although  more  slowly,  in  spite  of  the  presence  of  nonlinear 
elementary  divisors. 

2.k     The  Tchebychef  Iteration 

Because  of  the  recursive  property  of  the  Tchebychef  polynomials, 
the  third  criterion  from  Section  1.3  can  be  met.  An  iteration  based  on 
the  scaled  and  translated  Tchebychef  polynomials,  often  called  the  Stiefel 
iteration  (Stiefel  [19]),  can  be  carried  out  recursively,  requiring  the 
storage  of  only  three  vectors. 

As  before  let  x  be  the  solution  of  the  system  Ax  =  b.   Let  x 


n 


be  the  n^*1  iterative  solution,  let  e^  =  x  -x  be  the  error  at  the  n 

'      n      n 

step,  and  let  r  =  b  -Ax  be  the  residual  at  the  n^h  step.   If  the 

Tn(^) 
iteration  is  to  be  based  on  the  polynomials  P  (X)   =   ■= — ,  then 

n       T  (-) 
nvc 

e  =  P  (A)  e_  . 
n    nv  '      0 


We  have 


x   ,  -  x  =  (P  (A)  -  P  _ (A))  en 
ml    n     nv      n  i  1      0 


Let  Dx  =  x    -x  .   Since  the  Tchebychef  polynomials  satisfy; 


TQ(z)  =  1  , 


Tx(z)  =  z  , 


T  _(z)  =  2z  T  (z)  -  T  .  (z),  n  >  0  , 
n+lv        nv  '         n-lv  " 


then  for  n  >  0  we  have 


ViM  ■ 


C  n        G  n_1 


"n+1     c 

t   A-) 

n+lvc; 


"n-1     c 


T    £) 

nvc 


0\T(£±)        2  |  T    (^)        T     -  (^) 

2       n     c  c     n     c  n-1     c 

c 


Tx'(5)               T     A~) 
n+lvc'                 n+lvcr 

n+1 v  c ' 

T    (-)    " 
2       nV 

\  P   (X)    + 

n 

"2  *T  (*r 

c    n  c 

P  00  - 

n 

n-lxc 

T     ,  (-) 
.       n+lvcyJ 

T        (-) 
n+lV  _ 

T        (-) 
_    n+lV. 

p  noo 

n-1 


We  have  for  n  >  0 


P    00    -   P     . 00 

n     '  n+1 


"0     Tn(|) 

d       n  c 

X 

P  (x)  + 

n 

"       2  It  (*)  ' 

1            c    n  c 

T      ,£) 

_           n+1 v  c     _, 

P    00 
n 

C   T      .  (S) 

n+1 v  c ' 

+ 

"t    _(V 

n-1  cr 
^_  n+lvcy  _ 

p    n00 

n-1 

2       nV 

°  T        A 
n+lV_ 

X 

p  00  + 

n 

n-lV 

_Tn+iy_ 

(P 
n- 

a00  - 

30 


We  can  write 


Dx  =  (P  (A)  -  P  _  (A))  e„ 
n   v  nv      n+1 K    J  '      0 


"t  (*) 

nvc 

T  A±) 
n+l^cr 


A   Pn(A)  eQ  + 


T  _(*) 

n-1  c 

n+lV 


(P  .  (A)  -  P  (A))  e, 
v  n-1      nv  ''     ( 


Notice  that 


and  that 


r  =  A  e  =  AP  (A)  e_  , 
n     n     n     0  ' 


Dx  _  =  (P   '(A)-  -  P  (A))  e^ 
n-1   v  n-1       n      0 


We  have  then 

Dx  = 
n 

2 
c 

fT  (")  1 
n  c 

r  + 
n 

n-1  c 
n+1 v  c ' 

Dx  .  , 
n-1  ' 

T   ,A 
^  n+1  c  _ 

for  n  <  0.  For  n  =  0  we  have 


Since 


we  have 


DxQ  =  (PQ(A)  -  P1(A))  eQ  . 


PQ(\)  -  Px(\)  =  1  -•  (1  -  |)  =  \, 


n     1  a     1 

uxo  ~  I   o  "  d  ro  ' 


The  three  term  iteration  becomes:  given  initial  guess  x_  and 


parameters  d  and  c,  let 


ro  =  b  -  A  xo  ' 
Dxo  =  I  ro  > 


xi  -  xo  +  ^o  • 


For  n  >  0,  let 


r  =  b  -  A  x  , 
n         n  ' 

2   nV        n-lV  n 

Dx  = y-  r  +  3—  Dx  .  , 

n   c     (d}  n      A)       n-X 
n+1  c  n+1  c 


x  ,  =  x  +  Dx 

n+1    n     n 


The  coefficients  can  be  recursively  generated  and  the  iteration 

can  be  carried  out  with  only  x.,  r.,  Dx.  in  storage. 

i7   l    l 

If  the  spectrum  of  the  matrix  A  lies  in  the  right  half  plane, 
then  it  can  be  enclosed  in  an  ellipse  that  does  not  contain  the  origin  in 
its  interior.   The  associated  scaled  and  translated  Tchebychef  polynomials 
meet  the  criteria  established  in  Section  1.3.   They  have  minimal  maximum 
modulous  properties  on  ellipses,  their  derivative  sequences  also  converge, 
and  the  iteration  can  be  carried  out  by  a  three  term  recursion .   The 
remainder  of  this  thesis  will  be  devoted  to  implementing  an  iteration  based 
upon  the  Tchebychef  polynomials. 


3-   CHOOSING  OPTIMAL  PARAMETERS 

The  spectrum  of  the  matrix  A  can  be  enclosed  in  many  different 
ellipses.   In  fact,  given  any  family  of  ellipses,  ^(d,  c),  there  is 
some  member  of  the  family  that  contains  the  spectrum  of  A  in  its 
interior.   If  the  spectrum  of  A  lies  on  the  interior  of  F„(d, c),  the 
member  of  the  family  *^(d, c)  passing  through  the  origin,  then  the 
iteration  based  on  the  associated  scaled  and  translated  Tchebychef 
polynomials  will  converge.  We  would  like  to  choose  ^(d,  c)  so  that  this 
convergence  is  optimal  in  some  sense. 

In  Section  5-1  it  will  be  shown  that  the  choice  of  parameters 
which  yield  the  best  rate  of  convergence  can  be  found  as  the  solution 
of  a  mini -max  problem.   The  mini -max  problem  will  be  restricted  and 
restated  in  Section  5-2.   In  Section  3-3  it  will  be  shown  that  if  the 
matrix  is  real,  the  iteration  can  be  carried  out  in  real  arithmetic. 

3«1  The  Mini -max  Problem 

Suppose  d  and  c  have  been  chosen.   Let  F  (d, c)  be  the  smallest 
member  of  the  family  f(d,  c)  containing  the  spectrum  of  A  in  the  closure 
of  its  interior.   There  is  some  eigenvalue,  say  \   ,  such  that  \g  eF  (d, c) 
From  Section  2.1  we  know  that  F  (d, c)  is  associated  with  a  convergence 
factor,  r(\o),  and  that 

r(X   )  =  max  r(\  )  , 
S     >^i    X 


because  all  of  fcfae  other  >'i;  ftwalues  are  inside  or  on        (4»c) 
Recall  that  r(\)   is  also  a   function  of  d  and  c.     We  have 


r(M   -    |e 


[cosh"    (— — )  -cosh"    (-)  1 
v   c  c 


If  we  use  the  log  form  of  the  cosh"'  ,  then 


r(\)  = 


,d-\. 


,a-\, 


(§)  +  ((f)2-  l)5 


(d-\)  +  ((d-xf   -  c2)2 

1 
d  +   (d2  -  c2)2 


One  way  to  optimize  the  choice  of  d  and  c  is  to  make  r(\  )  as 
small  as  possible.   The  parameters  d  and  c  will  then  satisfy 


min  max  r(\.)  =  min  max 


d,  c  \^ 


d,  c  Xj_ 


1 
(d-\.)  +  ((d-*..)2  -  c2)2" 

d  +   (d2  -  c2)2 


A  more  rigorous  argument  which  yields  this  same  mini -max  problem 
is  as  follows.  With  a  polynomial  based  gradient  method  the  error  is 
suppressed  in  accordance  with  the  equation  e  =  P  (A)  e  .  The  following 
definition  of  rate  of  convergence  is  used  by  Young  [26]. 

Definition  The  rate  of  convergence  of  a  polynomial  based  gradient  method 
applied  to  the  system  Ax  =  b  is 


R(A)  =  -H   (lim(!|Pn(A)||n)) 


ii-k» 


7A 


We  would  like  to  choose  d  and  c  to  make  R(A)  as  large  as 

1 

possible  or,  equivalently,  to  make  lim  (|JP  (A)  ]|n)  as  small  as  possible 

n-»oo    n 

Let 

[cosh"  (-^-)  -  cosh"  (-)] 
M(\)  =  e         c  c   . 

If  we  use  the  log  form  of  cosh"  this  becomes 

1 

2    2N2 


m(M  =  ^  +  ^-^     -^  >  . 


d  -  c  )^ 


From  Lemma  2.2  we  have 

,n 


Pn(A.)  =  (M(\))J 

for  large  n.   Since  M(\)    is  analytic  in  an  open  set  containing  the 
spectrum  of  A,  there  exists  an  operator  M(A) .   The  eigenvalues  of  M(A) 
are  M(\j_)  where  \^   is  an  eigenvalue  of  A  (Dunford  and  Schwartz  [^]).  We 
have 

Pn(A)  =  (M(A))n 


for  large  n,  so  that 


1  1 

lim  (||P  (A)||n)  =  lim  (||M(A)n||n) 


n 
n->  co   ij-        n-*00 


=  spectral  radius  of  M(A) 
The  spectral  radius  of  M(A)  is 


max  |m(X.)|  =  max  r(\.) 


The  choice  of  d  and  c  which  yields  the  optimal  rate 

convergence  is  the  solution  to  the  mini -max  problem  above.  :;ince 

r(\)  is  a  function  of  d  and  c  as  well  as  \,  let  us  write  the  mini -max 

problem  as 

min  max  r(\. , d,  c)  . 
d,c  \± 

3-2  Restrictions 

In  this  section  we  will  show  that  if  A  is  a  real  valued 
matrix,  the  mini -max  problem  can  be  restricted  so  that  the  maximum  is 

taken  over  a  subset  of  the  eigenvalues  and  the  minimum  is  taken  over 

2 
d  and  c  such  that  d  and  c  are  real. 

In  Section  3.1  we  defined  the  ellipse  F  (d, c)  to  be  the 

s 

smallest  member  of  the  family  7"(cL,  c)  enclosing  the  spectrum  of  A. 
Since  Fs(d,  c)  is  convex,  the  eigenvalue  \  eF  (d,  c)  has  the  property 
that  it  is  on  the  boundary  of  the  convex  hull  of  the  spectrum.   In  fact, 
it  is  a  vertex  of  the  smallest  convex  polygon  enclosing  the  spectrum. 

Definition  Let  H  =  (\. |\.  is  a  vertex  of  the  smallest  convex  polygon 
enclosing  the  spectrum  of  A} .  We  will  refer  to  H  as  the  hull  of  the 
spectrum. 

The  elements  of  H  completely  determine  the  mini-max  problem. 

Lemma  3  •  1  For  any  d,  c 


max  r(X.,d,  c)  =  max  r(\.,d,  c) 
\±   '  X       \±e  H   '  X 


Proof  Suppose  there  is  a  d  and  c  and  X,    £H  such  that 

max  r(Xi,d,c)  =  r(^k,d,c)  . 


*-i 
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Since  X     /h,  we  know  that  X     is  not  a  vertex  of  the  smallest  convex 
polygon,  P,  enclosing  the  spectrum.   Let  F  (d,  c)  be  the  member  of 
0-(d,  c)  passing  through  X    .   Since  F  (d, c)  is  an  ellipse  and  passes 
through  a  point  that  is  either  in  the  interior  of  P  or  on  P  between 
vertices,  there  is  some  vertex  of  P  outside  F,  (d,  c)  .   From  the  definition 
of  P,  this  vertex  must  be  an  eigenvalue,  say  \..   By  Theorem  2.5 

J 

r(X  ,d,  c)  >  r(Xk,d,  c)  . 

This  contradiction  proves  the  lemma. 

If  A  is  a  real  valued  matrix,  then  the  eigenvalues  of  A  are 

real  or  appear  in  complex  conjugate  pairs.   The  hull,  H,  of  the  spectrum 

is  symmetric  with  respect  to  the  real  line.   This  motivates  the  following 

theorem,  the  proof  of  which  is  omitted. 

Theorem  j .2  If  A  is  a  real  valued  matrix  and  d  and  c  satisfy  the  mini- 

max  problem, 

min  max  r(>v..,d,  c)  , 
d,c  X±ett         X 

then  the  family  ^(d, c)  is  symmetric  with  respect  to  the  real  axis. 
Proof  Omitted. 

If  A  is  real  v;e  may  restrict  our  search  to  those  d  and  c 
which  correspond  to  families  of  ellipses  which  are  symmetric  with  respect 
to  the  real  line.   Such  a  family  has  foci  that  are  either  both  real  or 

are  a  complex  conjugate  pair.   Since  the  foci  are  d+c  and  d-c,  then  d 

2 
is  real  and  c  is  either  real  or  pure  imaginary.   In  either  case  c  is 

real.   Notice  in  the  log  form  of  the  definition  of  r(\,d,c)  in  Section  3.1, 

that  c  appears  only  as  c  .   Let  c2  =  c  ,  and  let 


r(\,d,c2)  = 


M-O  i    ((d-\)'  -  c2) 


d  +  (d  -  c2)' 


Since  the  families  ^(d,  c)  and  <^(d, -c)  have  the  same  foci,  d+c  and 
d-c,  the  parameters  d  and  c2  uniquely  determine  the  family  ^(d,  c) . 

For  A  real  we  may  restrict  d  and  c2  to  real  values.   In 
addition  we  may  ignore  those  values  of  d  and  c2  for  which  convergence 
clearly  does  not  occur. 

o 

Definition  Let  R  =  {(d,c2)/0  <  d,  c2  <  d^"}  . 

Corollary  3»3  If  A  is  a  real  valued  matrix  with  eigenvalues  in  the 

right  half  plane  the  mini-max  problem  can  be  written 

min   max  r(\.,d,  c2)  . 
(d,c2)eR  \.eH    X 

Proof  It  is  clear  from  the  discussion  above  that  d  and  c2  may  be 
restricted  to  the  real  numbers.   By  hypothesis,  A  has  eigenvalues  in 
the  right  half  plane.   If  d  <  0,  then  every  eigenvalue  of  A  would  be 

outside  Fn(d, c),  the  ellipse  passing  through  the  origin.   Convergence 

p 
could  not  occur  for  this  choice  of  d.   If  d  >  0  and  c2  >  d  ,  then  c  >  d 

and  d-c  <  0  <  d+c.   The  family  ^(d, c)  has  one  foci  on  each  side  of 

the  origin.   The  ellipse  F0(d, c)  is  the  degenerate  ellipse.   Since 

F0(d, c)  has  no  interior,  there  is  no  region  of  convergence. 

Because  A  has  its  eigenvalues  in  the  right  half  plane,  there 
is  some  d  and  c2  for  which  convergence  will  occur..   The  solution  of 
the  mini-max  problem  is  in  the  set  R,  and  the  lemma  is  proved. 

Notice  that  for  (d, c2)  e R  we  have 

r(\,d,c2)  =  r(\,d,c2)  . 
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The  maximum  is  completely  determined  by  the  eigenvalues  with 

nonnegative  imaginary  part. 

Definition  Let  H+  =  {A..eH/lm(A)  >  0]  .  We  will  refer  to  H+  as  the 

positive  hull  of  the  spectrum. 

Corollary  3 -^  If  A  is  a  real  valued  matrix  with  eigenvalues  in  the 

right  half  plane,  then  the  mini -max  problem  can  be  written 

min    max  r(A..,d,  c2)  . 
(d,  c2)eR  XieII+    1 

Proof  The  proof  is  clear  from  the  discussion  above. 

A  further  reduction  of  the  set  H+  is  possible.   We  would  like 
to  find  the  smallest  set  of  eigenvalues  that  completely  determines 


max  r(A..,d,  c2)  when  (d,  c2)  eR. 


A-i      x 

Definition  Let  K  =  [\.  ell+/there  exists  (d,  c2)eR  such  that  r(A..,d,  c2)  = 

max  r(\.  ,d,  c2)  J  .   The  elements  of  K  will  be  known  as  key  elements. 
X±  x 

Clearly,  if  (d,  c2)  eR,  then 


max  r(A..,d,  c2)  =  max  r(\.,d,  c2)  . 

X.eK    X  X.  €H^"    1 

l  i 

Criteria  to  determine  when  an  eigenvalue  is  in  the  set  K  are  needed. 
Lemma  j . 5  If  A.,  cK,  then  one  of  the  following  is  true: 

1.  Re  (A.,  )  <  Re  (A..)   for  every  \.   eH+ 

2.  Re  (A.,  )  >  Re  (A..)   for  every  A..  eH+ 

K   —       1  1 

3«   There  exist  A.^,  A.  eH+  such  that  there  is  an  ellipse, 

F,  (d,  c),  with  (d,  c2)  e  R,  passing  through  A.,,  A.  ,  and 

X  ,    containing  the  spectrum  of  A  in  the  closure  of 
its  interior. 


mot'     KV  .  (His  associated  with  a    family  o 

V    l,  c).     Let  Fw(d, c)   be  the  member  of   9^(d,c)  passing    through  >^. 

As    (d,  c2)    is  moved    through   the    rty.ion  R,    P'   (d,  c)    is   continuou 
deformed. 

If  X,    eK  then  there  is  some  point   (d  , c2   )  eR  such  that 

r(\k,d1,c21)   =  max  r(\i,d1,c21)    . 
Xi 

Since  r(\.,d  ,  c2  )  is  maximal  we  know  from  Theorem  2.5  that  F,  (cL,  c  ) 
contains  the  spectrum  in  the  closure  of  its  interior. 

Suppose  X     is  the  only  eigenvalue  on  the  ellipse  F,(d  ,  c,)  . 
If  X     does  not  satisfy  1.  or  2.  of  the  hypothesis,  then  there  are  at 
least  two  other  eigenvalues  in  H+,  say  X.    and  X   ,  such  that 

Re  (\  )  <  Re  (X^)   <   Re  (Xj  . 

Consider  deforming  the  ellipse  F^(d  ,  c  )  into  the  degenerate  ellipse, 
the  line  segment  connecting  X-^   and  X-^,    by  moving  (d,  c2)  through  R  (see 
Figure  3.1).   The  eigenvalues  X.   and  X     are  inside  Fk(d  ,  c,  )  but 
outside  the  degenerate  ellipse.   One  of  the  intermediate  ellipses 
must  have  passed  through  another  eigenvalue.  As  (d, c2)  moves  from 
(dn,c2,  )  let  fd0,c20)  be  the  first  point  such  that  the  ellipse,  F,  (d  ,  c  ) 
passes  through  another  eigenvalue,  say  X    .   Since  it  was  the  first, 
F,  (d  ,  c   )    still  encloses  the  spectrum,  and,  from  Corollary  ~5-h, 

Suppose  X,  and  X  are  the  only  eigenvalues  on  the  ellipse 
F^(d  ,  cp).  We  can  move  (d,  c2)  through  R  in  such  a  way  that  F-^(d,  c) 
passes  through  \y   anc^  *-p  •  As  c^  gets  negatively  large  the  foci  of  the 

\ 


1*0 


Xk   ^k^ci> 


i    /Fk(d2»C2) 


Figure  3-1 


Figure  3-2 


ellipse  F  (d,  c)  have  large  imaginary  part  and  the  ellipse  is  deformed 
into  the  infinite  column  between  Re(\)  and  Re(\»)  (see. Figure  3-2). 
If  Re(\,  )  <  Re(X.«)  then  one  of  the  intermediate  ellipses  must  have 
passed  through  \..      If  Re(\  )  >  ~Re(X   ),    then  one  of  the  intermediate 
ellipses  must  have  passed  through  X    .   In  either  case,  as  (d, c2)  moves 
from  (d^,  c2g)  let  F '  (d,,c;, )  be  the  first  ellipse  to  pass  through  a 
third  eigenvalue,  say  X    .   Since  it  is  the  first,  F,  (d-,,c^)  still 
encloses  the  spectrum,  and,  from  Corollary  $.k,   X     gH  .   This  proves 
the  lemma. 

The  lemma  provides  criteria  by  which  certain  elements  of  the 
set  H+  can  be  ignored.   The  implementation  of  these  criteria  will  arise 
naturally  from  the  algorithm  to  be  presented  at  the  end  of  Chapter  k. 

The  results  of  this  section  are  summed  up  in  the  following 
theorem. 


J.  I 


Theorem  5»6   IT  A  is  a  real  valuml  matrix  with  eigenvalues  in  the 

right  half  plane,  then  the  parameters  d  and  c  which  yield  the  optimal 

o 
rate  of  convergence  can  be  found  in  terms  of  d  and  c2  =  c  as  the 

solution  of  the  mini -max  problem 


min    max  r(\.,d,  c2) 
(d,c2)eR  \.eK 


3-3  Real  Arithmetic 

In  this  section  it  will  be  shown  that  if  d  and  c2  are  real, 
the  iteration  based  upon  the  associated  scaled  and  translated  Tchebychef 
polynomials  can  be  performed  in  real  arithmetic,  even  when  the  scaling 
parameter  c  is  pure  imaginary. 
Theorem  3-7  If  d  and  c2  are  real,  then 

T  (^) 

p.  00  = 


T  (*) 
n^c' 


is  a  polynomial  in  X   with  real  coefficients. 

Proof  Since  c2  =  c  ,  if  c2  >  0,  then  d  and  c  are  both  real,  and  7   (\) 

has  real  coefficients.   If  c2  <  0,  then  c  =  is  where  s  is  real.  We 

have 

T  £*) 

T  (t-) 
n  is 


If  n  is  even,  then  T  (z)  is  an  even  polynomial.  Each  term 

raises  (— r— )  to  an  even  power,  so  T„(— r— )   has  real  coefficients  and 
v  is           r    '     nv  is 

T„(t— )  is  real. 


k2 


If  n  is  odd,  then  T  (z)  is  an  odd  polynomial.  Each  term 


,d-\ 


raises  (— r— )  to  an  odd  power,  so  i  can  be  factored  out  of  each  term. 


is 


a  %  . 


Likewise,  T  (t— )  is  pure  imaginary.   The  quotient,  P  (M,  has  real 
coefficients,  and  the  theorem  is  proved. 

From  Section  2.5,  the  iteration  based  upon  the  scaled  and 
translated  Tchebychef  polynomials  is  as  follows:   to  solve  Ax  =  b, 
given  xQ,  d, c,  let 

rQ  =  b  -  AxQ 

Dx0  =  a  ro 
Xl  =  xo  +  Dxo  • 


For  n  >  1,  let 


r  =  b  -  Ax 


Dx  =  - 
n 


n+1 


2        nV 

°   T     _  (*) 
n+1  c  _ 

r     + 
n 

"t    A-)" 

n-1  cy 

T    A-) 

n+1   c  _ 

x     +  Dx 

n           n 

Dx 


n-1 


The  iteration  parameters  are  given  in  terms  of  d,  c  and  T  (— ) . 
If  c  is  pure  imaginary  the  computation  seems  to  require  complex 
arithmetic .   The  parameters  can  be  generated  recursively,  however, 
in  terms  of  d  and  c2. 

Theorem  3-8  If  d  and  c2  are  real  then  the  iteration  based  upon  the 
associated  scaled  and  translated  Tchebychef  polynomials  can  be  performed 
in  real  arithmetic. 
Proof  For  n  >  1  let 


Pl(n)  =  % 


T  (-) 
2   nV 


T  ,  (i)  ' 

n+1  c 


For   n  -  1  we  have 


T   A  * 

Pl(l)    =  -141  =  £  _£ £d_ 

T0(-)        C  2(-)    -1  2<T  -c2 

2  c  c 


T    (-) 
P2(D   =     °  (  1 


T2(§)        2(f)2-!        2d2-c2 


If  we  use  the  recursive  formulas   for  the  Tchebychef  polynomials    (see 
Section  2.1),    we  have   for  n  >  1 

Pl(n)   =  ^       n  c  c     nS 


c   _        ,dN        d     m        /dN 
T  ,,  (-)  T     ,  (-) 

n+lvc;  n+1  c 


1,(J)      ' 

n+lvc 


T     ,(-)  T     n(-) 

P2(n)   =  -iH  n"1  C 


=  -r(l  +  P2(n))    , 


T       '(*)       2  *  T  (*)   .  T       (1) 

n+lvcy  c     nvcy  n-1  c 


n-lvcy 

T        (-) 
2     n-lV 

T   $ 

nxc' 

C     T   (*) 

n  c 

T        (-) 
n  d         n-lV 

"  c     t  (i) 

n  cy 

4 

T       (-) 
d        2     n-lV 

<?    0  T  (a 

n  c 

PI (n-1) 

c2  PI (n-1) 

J5 

k  -^  -  PI (n-1)        4d  -  c2  PI (n-1) 
c^ 


induction 


kk 
If  d  and  c2  are  real,  then  Pl(l)  and  P2(l)  are  real.   By 


vo(„)  c2  Pl(n-l) 

P2(n)  =  kd  -   c2  Pl(n-l)  ' 


Pl(n)  =  |(1  +  P2(n)) 


are  real  for  n  >  1.   The  iteration  can  be  performed  in  real  arithmetic, 
which  proves  the  theorem. 


h.      SOLVING  THE  MINI-MAX  PROBLEM 

If  the  eigenvalues  in  the  positive  hull,  H+,  are  known,  then, 
as  was  shown  in  Chapter  3,  the  optimal  parameters  can  be  found  as 
the  point  that  minimizes  the  maximum  of  a  finite  number  of  real  valued 
functions  of  two  real  variables.  Consider  each  function,  r(A..,d,c2), 
to  be  a  surface  above  the  d, c2-plane.   Section  k.l   will  show  that  the 
mini -max  solution  is  either  a  local  minimum  of  one  of  the  surfaces,  a 
local  minimum  along  the  intersection  of  two  surfaces,  or  a  point  where 
three  surfaces  intersect.   Section  k .2  will  show  that  each  surface  has 
only  one  local  minimum  in  R.   It  will  be  found  explicitly.   The  local 
minimum  along  the  intersection  of  two  surfaces  will  be  found  in 
Section  k .J>   as  the  root  of  a  fifth  degree  polynomial.   In  Section  k.k 
the  unique  intersection  of  three  surfaces  will  be  found,  when  it  exists, 
and  existence  criteria  will  be  established.  An  algorithm  to  find  the 
mini -max  solution  among  the  possible  candidates  will  be  developed  in 
Section  4.5-   It  will  be  assumed  that  [\.)    are  eigenvalues  of  A  and 
that  Re(\.)  >  0,  ^(X.)  >  0. 

k.l     The  Alternative  Theorem 

The  following  theorem  will  enable  us  to  solve  the  mini -max 
problem  by  looking  at  three  or  fewer  functions  at  a  time. 
Theorem  k.l      (Alternative  Theorem)   If  (f.(x,y)}  is  a  finite  set  of 
real  valued  functions  of  two  real  variables,  each  of  which  is  continuous 
on  a  closed  and  bounded  region  S  and 


k6 
M(x,y)  ='  max  t±(x,y)    , 


then  M(x, y)  takes  on  a  minimum  at  some  point  (x  ,y  )  in  the  region  S. 
If  (xn, yn)  is  in  the  interior  of  S,  then  one  of  the  following  holds: 

1.  The  point  (x  ,y  )  is  a  local  minim-urn  of  f.  (x,y) 
for  some  i  such  that  M(xQ,y0)  =  f. (x  ,y  ). 

2.  The  point  (x  , yn)  is  a  local  minimum  among  the 
locus  ( (x,y)eS/f .  (x,y)  =  f.(x,y)}  for  some  i  and 
j  such  that  M(xQ,y0)  =  f^x^y^  =  f  (x0,yQ)  . 

3-   The  point  (x  ,  y  )  is  such  that  for  some  i,  j,  and 
.   k,  M(x0,yQ)  =  f.(x0,y0)  =  f.(x0,y0)  =  fk(x0,y0)  . 

Proof  Since  M(x,  y)  is  continuous  on  S,  it  takes  on  its  minimum  at  some 
point  (x  ,y„)  in  S.   In  the  following  all  neighborhoods  of  (x  ,y  ) 
will  he  understood  to  be  in  the  interior  of  S. 

If  there  is  a  neighborhood  of  (x„,yn)  such  that  M(x,y)  =  f . (x,y) 
in  the  neighborhood,  then  (x„, yn)  is  a  local  minimum  of  f. (x, y)  and 
M(x0,y0)  =  t^yj. 

Suppose  that  there  is  no  neighborhood  of  (x  , y  )  in  which 

the  maximum,  M(x,y),  is  determined  by  only  one  function  but  there  is  a 

neighborhood  in  which  M(x, y)  is  determined  by  two  functions,  say  f. (x, y) 

and  f.(x,y).   Since  f . (x,y)  and  f.(x,y)  are  continuous,  then  M(x  , Yq)   = 

£±(x0,y0)   =   f  (xQ,y0).   Clearly,  (xQ,y0)  is  a  local  minimum  in  the 

locus  t(x,y)eS/f .  (x,y)  =  f  (x,y)}. 
-*-         J 

Suppose  that  in  every  neighborhood  of  (xn,y^)  the  maximum, 
M(x,y),  is  determined  by  at  least  three  functions.   Since  they  are  all 
continuous,  there  are  three  functions,  say  f.(x,y),  f.(x,y),  and  f,  (x,y), 
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such  that  M(x0,y0)  =     ,yQ)  =  r.(xQ,y0)  .  fk<   ...  |.  Biifl  prov 
the  theorem. 

Since  the  spectrum  of  A  lies  in  the      half  plane,  it 
can  be  enclosed  in  an  ellipse  symmetric  with  respect  to  the  real  line 
that  does  not  include  the  origin.  Equivalently,  there  is  some  point 
in  the  region  R  such  that 

max  r(\.  ,  d,  c2)  <  1    . 

H     '  ' 

In  the  next  section  it  will  be  shown  that  for  each  \.,    r(\.,d, c2)  >  1 

i'    1       — 

on  the  boundary  of  R.   There  is,  therefore,  a  closed  and  bounded 
subregion  S  c  R  which  contains  the  mini -max  solution  in  its  interior, 
and  we  can  apply  the  Alternative  Theorem  to  the  mini-max  problem. 

k.2     Minimum  Point  of  a  Single  Function 

In  this  section  the  one  local  minimum  of  the  function 
r(\.,d,  c2)  in  the  region  R  will  be  found  in  terms  of  X. . 

Each  point  (d,  c2)  eR  is  associated  with  a  family  of  ellipses 
in  the  \-plane  whose  members  are  the  level  lines  of  r(\,  d,  c2).   Let 
F. (d,  c)  be  the  ellipse  in  this  family  passing  through  X.  .  From 
Theorem  2.5  we  know  that  r(\, d, c2)  takes  on  the  same  value  at  each 
\  eFj_(d, c).   In  particular,  consider  \     in  Figure  k.l.      Since 
(d-X  )2  >  c2  we  have 

(d-O  +((d->02-c2)2 
r(V.,d,c2)  =  r(\Q,d,c2)  =  — . 

d  +  (d2  -c2)2 
If  we  let  (d-\  )  =  a2,  then  every  point  \  =  x  +iy  on  the  ellipse 
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b^~~ 

^-Fi.(dt.c2.) 

Xo/ 

0 

d-c 

d 

d+c       y 

F. (d, c)  will  satisfy 


Figure  k . l 


a2 


y 


a2-c2 


Letting  X.  =  x.  +  iy.,  we  can  write 


r(\.,d,c2)  =  (a2)   +  (a2-c2) 


d  +  (d2-c2)2 


subject  to  the  constraint 


o        p 
,d-X.)      y. 

-  + 


a: 


a2-c2 


=  1 


(If  the  ellipse  Fi(d,c)  is  degenerate,  the  above  constraint  does  not 
hold.  If  y  ^  0,  the  degenerate  case  gives  a2  =  0.  If  y.  =  0,  the 
degenerate  case  gives  a2  =  c2.  In  any  case,  a2  varies  continuously 
with  d  and  c2.) 


is  clear  that  if  \.    La  in  the  right  hair  plane,  we 
pick  (d,c2)  ..  R  such  that  the  ellipse  F.  (d,c)  does  not  contain  the 
origin.   Equivalently,  there  is  some  (d, c2)  eR  such  that  r(\j,d,  c2)  <   1. 
The  following  lemma  will  show  that  we  need  only  look  in  the  interior  of 
R  for  the  local  minima  of  r(\.  ,d,  c2). 

Lemma  k.2     If  Re(\.)  >  0,  then  r(\^f6.,c2)   >  1  on  the  boundary  of  R. 
Proof  We  have  R=  ((d,c2)/0  <  d,  c2  <  d2}  .   From  Corollary  3.3  we  know 
that  r(\.,d, c2)  >  1  for  d  =  0  and  c2  =  d  .  From  the  constraint  equation 


we  have 


(d-x.)2  <  a2  , 

i   — 

y2  <  a2  -c2  . 


2 
If  c2  <  (d-x.)  ,  we  can  write 
—     l 


1 

2    ^2 


|  d-x  J  +  |((d-X,r  -  c2)^| 


1       -  v  i 


<  r(\.,d,c2) 


d  +  (d2-c2)2 


2         2 
If  (d-x.)  <  c2  <  d  ,    we  can  write 


| d-x.  |  +  y. 
1 1— -  <  r(\i,d,c2) 


d  +  (d2-(d-x.)2)2 


In  either  case  we  have 


lim   r(\.,d,c2)  >  1  , 
(d,c2)^oo    x 
(d,c2)eR 


which  proves  the  lemma. 
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Since  r(k.,&,  c2)  >  1  on  the  boundary  of  R,  which  would  cause 

the  iteration  to  diverge,  we  are  only  interested  in  local  minima  that 

occur  in  the  interior  of  R.   The  next  two  theorems  give  the  minimum 

of  r(\.,d,  c2)  in  terms  of  X.  .      Theorem  k.3   treats  the  case  X.    =   x.  +iy. , 

•  y.  4   0.   Theorem  UA  treats  the  case  X.    =   x.  . 
17  i  7  ii 

Theorem  ^ . $  If  A--  =  x.  +  iy.,  y.  4   °>  there  is  only  one  local  minimum 

of  the  function  r(A..,d,  c2)  in  the  region  R.   It  occurs  at  the  point 

2 
(d, c2)  =  (x.,-y.)_,  and  at  this  point 


2,        yi 


r(\  ,x  ,-y±)  =  - 


Proof  In  the  form 


r(X.  ,&,  c2) 


w  2   2n2 

Xi   (xi  +yi^ 


1         1 
>2)2  +  (a2-c2)2 


* 


1 

d  +  (d2-c2)2 


it  is  clear  that  r(X..,d,  c2)  is  a  continuous  function  of  d,  c2,  and  a2,  for 
values  of  d,  c2  in  R  and  continuously  differentiable  except  where 

a2  =  0, 
a2  -c2  =  0  . 

From  the  constraint  equation, 


(*-*i)     yf 


a2   +  a2-c2     ' 

it  is  clear  that  a2  is  a  continuously  differentiable  function  of  d  and 

c2  except  where 

a2  =  0, 

a2  -c2  =  0  . 


nee  yi  4   0  and  •>:  -  e2  >  yV,  we  have  a2  -  c2  >  0  in  R.  Nov.- 


d  =  x.   and  c2  <  -y .  . 

i  J  i 


This  is  the  ray  shown  in  Figure  k  .2  and  corresponds  to  the  degenerate 
ellipse  through  \. . 


C2 


C2  =  d' 


Ut.-yf) 


Figure  I4- .  2 


Any  local  minimum  of  r(X.,d,  c2)  must  occur  at  a  critical  point. 
A  critical  point  is  a  point  at  which  either  all  directional  derivatives 
are  zero,  or  some  directional  derivative  does  not  exist.  We  have  shown 
that  the  points  on  the  ray  above  are  critical  points.   It  will  be  shown 
that  these  are  the  only  critical  points  so  that  any  local  minimum  must  lie 
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on  the  ray.   It  will  then  be  shown  that  the  only  local  minimum  on 
this  ray  is  at  the  point  (x.,-y. ). 

First,  we  divide  R  into  two  regions:   d  >  x.  and  d  <  x. . 

Case  I:   d  >  x. . 

i 

2 
Let  c2  -  pd  .   If  p  <  1  this  describes  a  curve  through  R. 

Every  point  in  R  is  on  some  curve  c2  =  pd  ,  p  <  1.   It  will  be  shown 

that  there  is  only  one  possible  critical  point  on  this  curve. 

a)   If  p  >  0,  then  c2  >  0  along  the  curve  c2  =  pd2.   Let 

a2      1   d2 
6  =  — x.    y   =  —  =  -7T.  We  can  write 
H   c2'  '        p   c2 


1 

(P)2 

4 

(P- 

1 

■if 

1 

(yf 

+ 

(7- 

1 

■if 

The  directional  derivative  along  this  curve  becomes 


IP'          IP' 
2           12                1 

r' 

(P)2              (P-D2 
1                     1 

(yf  f  (7-D2 

r            p' 

2          11 

(P)2(P-D2 

For  r'  =  0  we  must  have  P'  =  0.   From  the  constraint  equation  we  have 


2(d-x.)    (d-x.)2       (y.)2 

P' =-s  P'  =  2pd  , 


3        (P)2       (P-D2 


so  that  p '  =  0  gives 


d-x. 


or 


a2  =  d(d-xi)  . 


Plugging  this  back  into  the  constraint  equation  we 

1 


c2  =  Pd2  , 


a2  =  d(d-x. ) 


This  is  the  only  possible  critical  point  along  the  curve  c2  =  pd% 
0  <  p  <  1. 

b)   If  p  <  0,  then  c2  <  0  along  the  curve  c2  =  pd^.  We  can 
write 


(■ 

1 

■vf 

+ 

(1- 

1 

■3)5 

(■ 

1 

+ 

(1- 

1 

■yf 

Taking  the  directional  derivative  along  this  curve  we  get 


r'  =  I 


1     1  ' 

(-p)2(l-3)2 


Again  we  have  that  r'  =  0  implies  3'  =  0,  and  again  the  only  point  on 

the  curve  satisfying  this  is 

2    2 
1   Xi  +  yi 

1 


5h 


c2  =  pd2  , 


a2  =  d(d-x.)  . 
v   1 


2 
c)   If  p  =  0,  then  c2  =  0  along  the  curve  c2  =  pd  .  We  can 


write 


r(\,,d,c2)  =  (a2 


i'  '      '  d 


Taking  the  directional  derivative  along  the  curve  we  get 


1  a2'   .  1   (a2) 


2     1   d     ,2 

2        d 
(a2)2 


The  constraint  equation  for  c2  =  0  becomes 


(d-x. )2  +  y2  =  a2  , 


so 


a2'  =  2 (d-x. ) 
l 


Substituting  this  in  the  equation  above  we  get 


r'  =  ^—  (d(d-x.)  -a2)  . 

(a2)2d2 


If  r'  =  0,  then  a2  =  d(d-x. ) .   Using  the  constraint  equation  for  c2  =  0, 
we  have 

d(d-x.)  =  (d-x.)2  +  y2  , 


or 


2 
x.  +  y. 

d  =  , 

x. 
l 


c2  =  0  , 


a2  =  d(d-x.)  , 


as  the  only  possible  critical  point  along  the  curve  c2  =  0. 

p 
On  each  curve  c2  =  pd  ,  p  <  1,  there  is  only  one  possible 

critical  point.   The  locus  of  all  of  these  possible  critical  points  is 

the  curve 

c2  =  d(d-ri)  , 


where 


and  on  this  curve 


2    2 
Xi  +yi 

71  =— x— 

l 


a2  =  d(d-x. ) 


Figure  4.3  shows  the  curve  c2  =  d(d-T))  .   The  possible  critical  point  on 
any  curve  c2  -  pd2  is  where  the  two  curves  intersect. 

Along  the  curve  c2  =  d(d-r)),  d  >  x.  we  can  write 


l 
1  •        1 

(d(d-X.))2  +  (d(d-X.)  -  dtd-T]))2 

r(X..,d,c2)  =  i 1—j 

d  +  (d2-d(d-T]))2 


1        1 

(d-x^   +  (r\-x±) 

1      1 

,2    ,  s2 
a  +  (ti) 
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c2  =  d(d-77) 


c2  =  pd' 


Figure  k.3 


The  directional  derivative  along  the  curve  is 


1     1 

1         tA2      2^ 
r(d  +n   ) 


r'    = 


(d-x.) 


1  1 

I^d-x.)2  +   (t^x.)2) 
,2 


1  1 

(d     +  t\   ) 


111 


1 


1 


1 


2        11 


1     1 


(d2(d%2)     -     (d-X.)2((d-X.)2+(T1-Xi)2)) 


(d2+T)2)2(d-X.)2d2 


Since  d  >  d-x.    and  r\  >  t\-x..  , 


then 


Ill      1     _      1 

d2(d       -  (d-x.)r((d-x.r  H  (r,-x. 


Wo  can  conclude  that  r'  >  0  for  d  >  x.  and  that  there  are  no  criti- 
points  in  the  subregion  d  >  x. . 

Case  IT :  d  <  x. . 
I 

p 
As  in  Case  I  we  look  along  the  curves  c2  =  pd  .   Taking 

directional  derivatives  along  the  curves  we  again  find  that  r'  =  0 

implies  a2  =  d(d-x. )  .   Since  a2  >  0  for  d  4-   x.  and  d  <  x.  in  this 
i  r     i         i 

subregion,  this  condition  is  never  met.  We  can  conclude  that  there  are 
no  critical  points  in  the  region  d  <  x . . 

We  have  now  isolated  the  possible  local  minima  to  the  line 

d  =  x. •   From  the  constraint  equation  we  see  that  if  d  =  x. ,  then 
i  i' 


a2  = 


2  2  2 

y.  +  c2   for   x.  >  c2  >  -y.  , 

i  l         i 


2 
for  -y.  !  >  c2  . 


2  2 

Case  I:  x.  >  c2  >  -y. ,  d  =  x. . 


We  can  write 


1 
(yi  +  c2)   +  yi 

I 
2    ^2 


r(\,,d,c2)  - 


x. 


+  (x.  -  c2)' 


Taking  the  derivative  with  respect  to  c2  we  get 


1 


| 1   (x±  +(X2  -c2)2)  +  |  1 ((y^  +C2)2  +y. 


l         i' 


(y?  +c2)2                   (x2  -c2)2 
r'  =  = 1 


(x.  +  (x2  -  c2)2)2 
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Since  each  term  is  positive,  then  r'  >  0,  and  r(\.,d,  c2)  is  increasing 

2 
as  c2  increases  from  c2  =  -y.  along  the  line  d  =  x. . 

i  l 

2 
Case  II:   c2  <  -v.,  d  =  x. . 
i'  l 


We  can  write 


1 
r(\,,d,<£)  = 


1'   '  _L 

x.  +  (x.  -c2)2 

1   v  1 


Taking  the  derivative  with  respect  to  c2  we  have 


1  1 

2  x  (*±  +  (xi  ~c2)  )  "  g  1  (~c2) 

(-c2)2  (x2-c2)2 

r'  =  1 

1 

(x.  +  (x2  -  c2)2)2 
l     l 


1  1 

2      2         2      £" 
(x.  -c2)  (x.  +  (x7  -c2)  )  +  c2 

I  I     I 

2(x.  +  (x2  -c2)2)2(x?  -c2)2(-c2)2 
l     i  l 


1 

/  2   o\2    2 
x.  (x.  -  c2)   +  x. 

li     l 

1  1     1 

2(x.  +  (x?  -c2)2)2(x.  -c2)2(-c2)2 


x. 

i 


1        11 
2(x.  +  (x2-c2)2)(x.-c2)2(-c2)2 


nee  each  beno  is  positive,  then  r'  <  ,  and      L,c 


as  c2  increases  toward  c2  =  -y.  alotv,  t.he  line  d  =  • 

We  can  conclude  that  the  only  local,  minimum  of  r(\,  . 

2 
occurs  at  the  point  (d, c2)  =  (x.,-y.)  and  at  this  point 


r(\i,xi,-yp 


xi  +  ^x'i  +  yi^2 


This  proves  the  theorem. 

Theorem  h .4  If  \.  =  x.,  then  there  is  only  one  local  minimum  of  the 
function  r(\.,d,  c2)  in  the  region  R.   It  is  at  the  point  (d, c2)  =  (x.,0), 
and  at  this  point 


r(\,,x,,0)  =  0  . 


i'  i 


Proof  Since  X-    =  x. ,    we  have  a2  =  (d-x.)   except  for  the  degenerate 

2 
case.   The  degenerate  case  occurs  when  c2  >  (d-x.)  ,  so  that  X.  =  x.  is 

_  x       '  11 

2 
between  the   foci,    d+e   and  d-c.      If  c2  ;>   (d-x.)    ,    then  a2  =   c2.     We  can 

write 


^   |d-x.  |    +    ((d-x.)2  -c2)2 


r(\.,d,c2) 


d  +   (d2  -c2)2 


- 


(c2) 


P  ? 

d  +    (d~  -c2) 


for      (d-x. )      >  c2  , 

l  ' 


for      (d-x. )      <  c2    . 
i'      — 


In  this   form  it  is   clear  that   r(\.,d, c2)    is  continuous   in  R 
and  continuously  differentiable  except  where 


6o 


c2  =  "(d-x.)   , 


and  where 


d  =  x.   for  c2  <  0 

1 


As  in  the  proof  of  Theorem  k.3,    it  will  be  shown  that  the 
only  critical  points  of  the  function  r(\.,d, c2)  in  R  are  those  described 

above.   Divide  R  into  subregions  as  shown  in  Figure  k.k.      Region  I  is 

2 
where  c2  >  (d-x.)  ,  which  corresponds  to  the  degenerate  ellipse.   In 

2 
this  region  we  have  a2  =  c2.  Region  II  is  where  c2  <  (d-x.)  ,  d  >  x., 

2 
and  Region  III  is  where  c2  <  (d-x.)  ,    d  <  x. .   In  these  regions  we  have 

a2  =  (d-x.)2. 
v   l 


C2  =  d2    c2  =  (d-xL)J 


Figure  k.k 


Case  l:   c2  >  (d-x.)2. 


In  this  subrty.ion  we  have 


r(\«,d,c2)  = 


d  +  (d2  -c2)2 


Taking  the  partial  derivative  with  respect  to  c2,  we  get 

1  1 

\  -i-  (d  +  (d2  -c2)2)  +  \  1 -  (c2)2 

.,      (c2)" (d"-c2)" 

1 

(d  +  (d2-c2)2)2 


1         1 

1  (d  +  (d2  -c2)2)(d2  -c2)2  +  c2 

2  III 
(d  +  (d2  -c2)2)2(c2)2(d2  -c2)2 


Since  each  of  the  terms  is  positive,  we  have  r'  >  0  in  this  region  which 
implies  that  there  are  no  critical  points  in  this  region. 


Case  II:   (d-x.)  >  c2,  d  >  x. 
i  l 


Case  III:   (d-x. )^  >  c2,  d  <  x. . 
v   x      '      i 

These  two  cases  will  be  treated  together.   In  these  regions  we 
have  a2  =  (d-x.)  .  As  in  the  proof  of  Theorem  k.3,    consider  the  curves 
c2  =  pd2,  p  <  1.   Let 
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By  the  same  argument  used  in  Theorem  h.3,    the  directional  derivative 

2 
along  the  curve  c2  =  pd  is  zero  only  if  P'  =  0.  We  have 


P'  = 


2x. (d-x.) 
l    l 


pd 


3 


Clearly,  p1  /  0  in  either  subregion  II  or  subregion  III,  so  we  may 
conclude  that  there  are  no  critical  points  in  subregion  II  or  subregion  III, 
We  have  isolated  the  possible  local  minima  to  the  curve 


c2  =  (d-x.)   , 
i 


and  the  ray 


d  =  x.,  c2  <  0  . 
l 


Along  the  curve  c2  =  (d-x.)^  we  have 


r(\,,d,c2)  = 


d-x. 
1    l 


d  +  (d2  -  (d-x. )2)2 


If  we  take  the  directional  derivative  along  this  curve,  we  find 


-  (d  +  (d2  - (d-x. )2)2)  +  (d-x.)(l  + 


X. 

r> 


(d2-(d-x.)2)2 


for  d  <  x., 


- 


(d  +  (d2-(d-x.)2)2)2 


1 
,2^2 


(d  +  (d^  -(d-x.)2)2  -  (d-Xi)(l  +  1 1) 


(d2-(d-x.)2)2 


v. 


1 

(d  +  (d2-  (d-Xi)2)2)2 


for  d  >  x., 
l 


f 


-x. 


2nSWj2 


,2x2 


- 


(d   i    (d-  .(tij)T)(4'.((Ui)') 


X. 

1 


V 


1  1 

>x2Wj2      ,_,        ,2N2 


(d  +   (d2  -(d-x.)2)2)(d2  -(d-x.)2); 


for  d  <  x., 


for  d  >  x. 


Since  all  of  the  terms   are  positive  we  have 


and 


r'   <  0,      for     d  <  x., 

'  i' 


r'    >  0,      for     d  >  x.  . 

'  l 


The  only  local  minimum  along  the   curve   c2  =    (d-x. )      occurs  at  the  point 
(d,c2)   =    (x.,0). 

On  the  ray  d  =  x . ,    c2  <  0  we  have 


r(\i,d,c2) 


:-c2)' 


i 

2  2 

x.    +    (x.   -c2) 
li 


Taking  the  directional  derivative  along  the  ray  we  have 


1         1 

2   M2,    1,   „N2    -1 


l-^y  (x.  +  (x^-c2)")  -  |(-c2) 


(-C2)2 


r'  = 


(x2-c2)2 


1 
(x.  +  (x2-c2)2)2 


X. 

1 


1  1       1   ' 

(x.  +  (x2-c2)2)(x2_c2)2(.C2)2 


6k 


Clearly,  we  have  r'  <  0  along  the  ray,  so  that  r(\.,d,  c2)  decreases 
as  c2  increases  to  c2  =  0. 

We  can  conclude  that  the  point  (d, c2)  =  (x.,0)  is  the  only 
local  minimum  of  the  function  r(\.,d, c2)  in  the  region  R.  At  this 
point  we  have  r(\.,x.,0)  =  0,  and  the  theorem  is  proved. 


k.~5     Pair-wise  Best  Point 

In  this  section  it  will  he  shown  that  the  locus  of  points 
for  which  two  surfaces  intersect  is  a  continuous  curve  through  the 
d, c2-plane .  In  general,  the  minimum  along  this  curve  can  be  found  in 
terms  of  the  root  of  a  fifth  degree  polynomial  in  a  certain  interval. 
The  coefficients  of  this  polynomial  will  be  given.  It  can  also  be  shown 
that  this  minimum  is  the  only  local  minimum  along  the  curve,  but  it  will 
not  be  proved  here.  A  special  case  will  be  treated  in  which  the  minimum 
can  be  found  in  terms  of  the  root  of  a  third  degree  polynomial. 

Consider  all  points  (d,  c2)  €R  such  that 


r(\  ,d,c2)  =  r(\ ,,d,c2)  . 

From  Theorem  2.5  we  know  that  \.    and  \.   are  on  the  same  member  of  the 

i      3 

family  <p(d, c).   Moreover,  they  both  satisfy  the  equation  of  the  ellipse. 

If  X.  =  x.  +iy.  and  \.   =  x.  +iy.,  then  we  can  write 
iij      3  3  3' 

2     2 
(d-x  r    77 

+  ^"TT  =  1 


a2     a2-c2  ""   ' 
and 

(d-x.)     y 

a —  +  a—  =  1  . 

a2     a2-c2 


-pose  y.  ■  y.j  then,  we  have 


which  gives 


Let 


(d-x.f  =  (.i-x.r   , 


X  .  +x. 

d  =  -^— L 


X.  -X. 

A  _  J I 


X  .  +X. 


yj  'yi 

S  =  -jL 


2    ' 


m   YJ+yi 


2 


We  may  assume  that  x.  >  x.  ;  then,  we  have  A  >  0,  B  >  0,  T  >  0.   For 
S  =  0  we  have 


d  =  B  , 
and  the  ellipse  equation  becomes 


a2  +  a2-c2 


Solving   for  c2,    we  get 


a2(a2  -  (A2  +T2)) 
c2  =  — s s — '— 

(a2  -A^) 


Suppose  that  y.  4-  Y-j    that  is,  S  /  0.   Let  p  =  —  .   The  ellipse  equations 

1      J  Cc_ 

become 
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a         \2  2 

d_Xi}  yi  o 


and 

(d-x.)  y. 


+  tA-  =   c2    . 


P  3-1 


Solving   for  £,    we  get 


(x.-x.  )(2d-(x.-x. )) 
P=   (x.-x.,(2d-(x.-Xl,,   +  (y.-y.)(y.+y.) 


In  terms  of  A,    B,    S;    and  T,    this  is 

(d-B) 


(d-(B  +^)) 


Substituting  this  back  into  the  ellipse  equation,    we  get 


(A  \2  2 

(d-x   )  y 

c2  = —  +  — - 

P  0-1 


-(d-(B-A))2         |  (T-S)2 

(d-B)  |S 


S^ 


<a-(B  +  T»  (d-(B+f)) 


ST 

(d"(B+~})    ffd      fBAxx2    ,    A(d-B)(T-S)2 
jj-^  ((d-(B-A))      +  — ) 


(d-(B+f^))(d-(B   -   A|))(d-(B  -  4)} 
(d^Bl 


Notice  that 


so  that 


c2      l-(B  -  A|))(d-(B  -  A§)) 


=  §   (d-(B  -  A|))(d-(B  -  A§))  , 


a2  =  (d-(B  -  A|))(d-(B  -  A|))  . 


For  S  /  0  we  can  express  c2  and  a2  in  terms  of  d.  For  S  =  0 
we  have  d  =  B  and  can  express  c2  in  terms  of  a2.  To  insure  that  these 
parameters  describe  an  ellipse,  we  must  have 

a2  >  0  , 
and 

a2  -c2  >  0  . 

Lemma  k- . 5  The  above  equations  for  c2  and  a2  describe  the  parameters  of 
an  ellipse  passing  through  \.    and  \.   if 

d  <  B   when   S  <  0  , 
d  >  B   when   S  >  0  , 
and 

d  =  B,  a2  >  A2   when   S  =  0  . 

Proof  We  have  seen  that  d  =  B  when  3=0.   For  S  =  0  we  have 

a2(a2  -  (A2  +T?)) 
(a2  -A2) 


2 
a2  -A 
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2  2 

For  a2  >  A  we  have  a2  -  c2  >  0.  Notice  that  as  a2  approaches  A,  c2 

approaches  -co  and  when  a2  gets  large  positively,  c2  also  gets  large 
positively;  thus,  c2  takes  on  all  possible  values. 
When  S  /  0  we  have 

(a-(B  +S) 

c2  ■ a£  — r^Bi — ' 


so  that 


(d-(B  +S)) 
a2-c2  =  a2  -  a2  ^ 


ST 
A 


=  a2TaiBy 


If  S  <  0  and  a2  >  0,  then  we  must  have  (d-B)  <  0.   If  S  >  0  and  a2  >  0, 
then  we  must  have  (d-B)  >  0. 
For  a2  we  have 

a2  =  (d-(B  -  A|))(d-(B  -  A|))  . 

Since  A,  B,  T,  >  0,  then  if  S  >  0  and  d  >  B,  we  have 

(d-B)  +  a|  >  0  , 
(d-B)  +  a|  >  0  , 

so  a2  >  0.   If  S  <  0  and  d  <  B,  then  we  have 

(d-B)  +  a|  <  0  , 

(d-B)  +  a|  <  0  , 
so  a2  >  0.   This  proves  the  lemma. 


A  geometric  interpretation  of  this  result 
Figure  h.5.      If  S  >  0,  dun  y  .  ■  .v.  and  any  ellipse  throuj      <nd 
\.  must  have  its  center,  d,  closer  to  x.  than  x.,  or  B  <  d.   A 
similar  explanation  hoLds  for  S  <  0  and  S  =  0.   Figure  tiat 

the  locus  of  points  in  the  d,  c2-plane  might  look  like  in  each  of  t.he 
three  cases. 

Our  object  is  to  find  the  minimum  of  r(\.,d, c2)  =  r(\.,d,  c2) 
along  the  curves  shown  in  Figure  k. 6.  Theorem  k.6  will  treat  the  case 
S  3/  0,  and  Theorem  k.J   will  treat  the  case  S  =  0. 

The  following  notation  will  be  convenient.   Let 

z  *  d  -B; 
then,  we  can  write 

(z  +  A|)(z  +  A|)(z  -  21) 


c2 


a2  =  (z  +  A|)(z  +  A|)  . 


The  parameters  d,  c2,  and  a2  describe  an  ellipse  through  X.    and  \.   if 

z  >  0   for   S  >  0  , 
z  <  0   for   S  <  0  . 

Theorem  k.6  If  S  £  0,  the  point  at  which  r(\.,d, c2)  is  minimal  along 
the  locus  of  points  for  which  r(\.,d,  c2)  =  r(X..,d,  c2)  can  be  found  as 
the  root  of  the  polynomial 

5   ,  _  h-   .  _  3   ,  _  2 


plz  +  P2Z  +  P3Z  +  P4Z  +  P5Z  +  p6 
in  the  interval 


=  0 
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iy 

ii 


iy 

i 


-»»-*-x 


B  d>B 

S>0 


d<B 


B 
S<0 


*l 


B        d  =  B 
S  =  0 


-  x 


Figure  1+.5 


c2 


c2  =  d' 


S  >0 


S  <  0 


s  =  o 


Figure  k.6 
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.A)  S  > 

(-A,  for 


The   coefficients   are 


Pl  =    (2B  -  A(|  +  |))(2B  +^-  A(|h   |))    , 


p2  =    (2E  +  ^  -  A(|  +  |))((2AB  +ST)(|  4  |)    +  1+A2) 
+   B2(2B  -  A(|  +  |))    4    B(B2-A2)    , 


P3   =   hAk   -   WB(|  +  |)    +  A2GT((^  +  ^ )    -   5(|  +  |)) 

2         ,2 
f  A2B2(^  +  %  +  3)    , 
S  T 


p^   =  AST((B   -   A|)(B   -   3A|)    +    (B   -   A|)(B   -   3a|)  )    , 


p5   =   -3A'ST(2B   -  A(|  +  |))    , 


p6   =    -3A3ST(B2  -A2) 


Proof     We  will  refer  to  r(\.,d, c2)    as  r(z).     We  can  write 


i  (z   -   SZ)      i 

((z   +  A|)(z   +  A|))2   4-    ((z   +  A|)(z   +  A|)(l   - ^-))2 

r(z)    =  , . _ 

?        (z  +A|)(z   +a|)(z   -  ^)    2 
(z+B)    +    ((z+B)2 2 -1 A_) 


Since  A,  B,  T  >  0,  then  r(z)  is  a  continuously  differentiable  function 
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of  z  if  z  >  0  when  S  >  0  and  z  <-0  when  S  <  0.  Let  p 
for  c2  /  0.  We  can  write 


a2 
c2 


>    7  = 


c2 


r 


r(z)   =    < 


1             1 

(p)2  +  (p-i)2 
1               1 

for   c2  >  0  , 

(yf  +  (7-D2 

i              l 

(-P)2  +  (i-p)2 
l              i 

for  c2  <  0  , 

(-7)2   +    (1-7)2 


^ 


with  a  removable   singularity  at  c2  =  0.      Taking  the  derivative  with 
respect  to   z  we  have    for  c2  >  0 


r'(z)   = 


1  ft1        1       p] 

2  1  +2  1 

(P)2  (P-D2 


1  1 

((7)2+(7-D2)' 


1  1 

((yf  +  (7-D2)2 


2  12  1 

2  2 

(7)  (7-1) 


1  1 

((3)2+(3-l)2) 


1 


1  (P)     +  (P-D 

2  1  1 

(yf  "4  (7-D2 


P' 


(pro-D-      (7)  (7-1)' 


r(z)  /  p' 


(P)    (P-D  (7)    (7-D' 


and  for  c2  <  0  we  have 


r'(z) 


r(z) 


-P' 


(-P)"(l-P)"        (_7)"(i-7)' 


with  q  removable  singula  In  1  we  hn 


•B)2 
7  = 


(z  +A§)(z   4   A§)(z   -  fO 


ST 
fl  1    -  A 

~stT    » 


(z+B)2z   -    (z   +  A|)(z   f  A|)(z   -  |£) 
"  (z   +  a|)(z   +  A|)(z   -  §0 

Q(z)     


(z   +  a|)(z   4   a|)(z   -  |^) 


where  Q(z)  is  a  quadratic.   Taking  derivatives,  we  have 


P' 


ST 
A 


,     STX2   » 


((z+B)2  +2(z+b)z)(z  +  A|)(z  +  a|)(z  -  ^) 


7  = 


(z  +A|)(z  +  A§)(z  +f) 


O 


STs   ,     «SW     ST> 


-  z(z+B)^((z  +  A±)(z  4  A£)  4  (z  4  A±)(z  -  ^L)  +(z  4  A§)(z  -  |i)) 


C(z)(z+B) 


(z  +a|)(z  4  A|)(z  -  ^) 


7^ 


where  c(z)  is  a  cubic.  From  here  on  we  will  assume  S  >  0.   The 
arguments  are  identical  for  S  <  0.   Plugging  the  above  expressions 
into  the  equation  for  r'(z),  we  get 


r  (z)  =  - 


/_  i 


2 


(^)5  + 


A  '  1 


(z  -  f)(zf    \       (Q(z))2(z  +  a|)(z  +  A|) 


ST 
This  expression  has  a  singularity  at  z  =  j  which  corresponds  to  c2  =  0 

To  see  that  this  is  a  removable  singularity,  notice  that 
^,STN     /STWt3   STwST   ATwST   ASx 

n/STN    ST,ra   STx2 
Q(-)  =  ~(B  +  T)   , 


so  that 


(m f  + r_l^ .  o 


We  want  to  show  that  r'(z)  has  a  root  in  the  interval  (0,A) 
At  z  =  0  we  have 

C(0)  =  -ABST, 
Q(0)  =  AST  . 

Taking  the  limit,  we  have 


At    ■       A  we  have 


p  p 

C(A)    -    ,V   "/;'   '     (2(A2-ST)    -   A(AiB))    , 


Q(A)    =  A   (ST(A+B)2   -    (T+S)2(A2-ST))    , 


so  that 


(A2-ST)     P'  i 

(Q(A))^A2(S+T)^ 

I       /       i 

r(A)        (A;  ffSTs2         (2  (A   -ST)    -   A(A+B)) 

2        (A2-ST)     TA)       h  | 

\  (Q(a))2 


r(A)      (ST)"     /1  (2(A2-ST)    -  A(A+B)] 


Op  i 

(A   -ST)    \  g     p  i 

(ST(A+B)  "   -    (S+T)    (A~-ST)) 


We  would  like  to   show  that  r'(A)   >  0.      Removing  the  terms  known 
to  be  positive,    we  would  like   to   show  that 


^1 L    4  (2(A2-ST)    -  A(A+B)) 

(A2-ST)    I  \ 

(ST(A+B)^   -    (S+T)2(A":-ST))^ 


which  is   equivalent  to 

1 
A(A-fB)    -   2(A2-ST)         (ST(A+B)2   -    (S+T)2(A2-ST) )2 
(A2-ST)  (A2-ST) 
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If  (A  -ST)  >  0,  this  is  equivalent  to 


A(A+B)  -  2(A2-ST)  >  (ST(A+B)2  -  (S+T)2(A2-ST) )2 


If  we  square  both  sides  we  have  the  equivalent  expression 


A2(A+B)2  -  U(A+B)(A2-ST)  +  l4-(A2-ST)2  >  ST(A+B)2  -  (S+T)2(A2-ST)  , 


which  is  equivalent  to 


(A  -ST)(^(A2-ST)  -  U(A+B)  +  (B+k)d   +  (S+T)2)  >  0  , 


which  is  equivalent  to 


(A2-ST)((B-A)2   +    (T-S)2)    >0     , 


which  is  clearly  true. 

p 
If  (A  -ST)  <  0,  we  must  show 

1 

A(A+B)  -  2(A2-ST)  <  (ST(A+B)2  -  (S+T)2 (A2-ST) f      . 

Since  A(A+B)  -2(A2-ST)  >  0  when  (A^-ST)  <  0,  then  squaring  both  sides 
and  rearranging  as  before  gives  the  equivalent  expression 

(A2-ST)((B-A)2  +  (T-S)2)  <0  , 


which  is  clearly  true. 

Since  r'(0)  =  -co,  and  r'(A)  >  0,  and  r(z)  is  continuously 
differentiable,  we  can  conclude  that  r(z)  has  a  minimum  in  the  interval 
(0,A).  At  this  point  we  must  have  r'(z)  =  0.   From  the  expression  for 
r'(z)  we  must  have 
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or 


" rSM =  o 

(QU))2(z  +  a|)(z  +  a| 


(^5(Q(Z))5  = 


> 


A  (z+a|)(z+a|) 


Squaring  both  sides,  we  must  have 


^Q(z)(z  +  A|)(z  +  A§)  =  C(z)2 


Since  C(z)  is  a  cubic  and  Q(z)  is  a  quadratic,  this  gives  rise  to  a 

sixth  degree  polynomial.   One  of  the  roots  of  this  polynomial  is  at 

ST 
z  =  -jr  i    which  corresponds  to  c2  =  0,  an  extraneous  root.   Factoring 

out  this  root,  we  get  the  fifth  degree  polynomial  whose  coefficients 

were  given  in  the  hypothesis.   This  polynomial  must  have  a  root  at  the 

point  z  e(C,A)  at  which  r'(z)  =  0.   For  S  <  0  the  interval  is  (-A, 0). 

'This  completes  the  proof  of  the  theorem. 

It  can  be  shown  with  arguments  based  on  the  derivatives  of 

the  fifth  degree  polynomial  that  this  is  the  only  Local  minimum  of  r(z) 

We  next  turn  our  attention  to  the  case  S  =  0.   The  following 

notation  will  be  convenient.   Let 

y  =  a2; 
then,  we  can  write 

d  =  B  . 


r<?     y(y-(A2+T2)) 
(y-A2) 
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yT2 
a2  - c2  =  * 


(y-A2) 


The  parameters  d,  c2,  and  a2  describe  an  ellipse  through  the  points 

X.    and  \.   if  y  >  A2. 
■  l      J 

Theorem  h.^  If  S  =  0,  the  point  at  which  r(\.,d, c2)  is  minimal  along 
the  locus  of  points  for  which  r(\.,d,  c2)  =  r(X.,d,  c2)  can  be  found  as 
the  root  of  the  polynomial 

qxr  +  ^y  +  q^y  +  q^  =  0 
2  2 

in  the   interval    (A  ,  B   ).      The  coefficients   are: 


qx  -    (B2+T2)    , 


2   2 

4   2 
q^   =   3AV   , 

q^   =   -AifB2(A2  +T2) 


Proof  We  will  refer  to  r(\.,d,c2)  as  r(y) .  We  can  write 

r(y)  . iz=L 


1 

lB2.y(y-(A2+T2))\ 
(y-A2) 


2  2 

Notice  that  for  y  e  (A  ,B  ),  r(y)  is  a  continuously  differentiable 

a2 

function  of  y.   If,  as  in  the  proof  of  Theorem  h.G,    we  let  p  =  — ^  , 

d2        / 
7  =  —  for  c2  f  0,   we  have 


r 


r(y) 


] 

1         T  "        \_         1 

(7)2(7-D:' 


, 


r'(y)  ■   < 


r(y) 


11  11 

(-^)2(i-3)2       i-yfd-y)2 


for  c2  >  0   , 


with  a  removable   singularity  at   c2  =  0.      In  terms  of  y  we  have 


3  = 


(y-A2) 


(y-(A2+T2))  ' 


B2(y-A2) 


y(y-(A2+T2)) 


3-1   = 


(y-(A2^T2)) 


a  =   -(y2  -(A2+B2+T2)y-fA2B2) 
y(y-(A2  +T2)) 


Taking  the  derivative  we  have 


P' 


■r2 


(y-(A^4-T2))2 


7'    = 


B2(y2  -2A2y  +A2(A2  +T2)) 

(y(y-(A2+T2)))2 


8o 
Combining  the  equations  above,    we  have 


r  .  (y)   =  iM  I (        B(y2  -2A2y+A2(A2  +  T2))  _   , 

(y  -  (A2  +T2))(y-A2)2    \y(-(y2  -  (A2+B2+T2)y  +A2B2))2 


We  would  like  to  show  that  r'(y)  has  a  root  in  the  interval 

2  2 
(A  , B  ) .  We  can  write 


t    M    \        r(A^)    i-     1      i    ^B-A^^ 
lim  r'(y)  =  -A^_  .  lim .  (.(_)) 

y-^A^  y-A^      « 

(y-A2)2 


If  we  let  y  =  B  ,  we  have 


r.(B2)  =  £i£l 1 /b((b2-a2)2+a2t2  _  T 


(B2  - (A2  +T2))(B2-A2)2  \   B2(B2T2)2 


r(B2) 1 /  (B2-A2)(B2  -(A2  +T2)) 

2 
B  T 

(B2  -(A2  +T2))(B2-A2)2  ' 


1 
r(Bg)  (B2-A2)2 
2     B2T 


Since  lim  r'(y)  =  -«,  r'(B  )  >  0,  and  r'(y)  is  continuously 

y-A2 

2  2 
differentiable  in  the  interval  (A  ,B  ),  then  r(y)  must  have  a  minimum  in 


2 
the  interval    (A' ,  B'  K      From  the  equation  '(y)   we  see  thm 

re  have 


-2A2y  +AP(A2  4 T2))  =  T  y(-(y2  -  (A2+B2+T2)y  +A2B2))2 


If  we  square  both  sides  of  this  equation  and  rearrange  terms,  we  have 

o   p 
a  fourth  degree  polynomial  in  y.   This  polynomial  has  a  root  at  y  =  A  +T  , 

which  corresponds  to  c2  =  0.   Factoring  out  this  extraneous  root  we  have 

the  third  degree  polynomial  whose  coefficients  were  given  in  the  hypothesis 

Any  point  y  such  that  r'(y)  =  0  will  be  a  root  of  this 

polynomial.   If 


then 


3  2 

Q(y)  -  qxy    +  q^y    +  q^y  +  q^  , 


Q1  (y)  =  3qxy    +  ^y  +  q5 


Taking  the  discriminant  of  Q'(y),    we  have 


P  h.   h.  li   P      P        P 

l*q|    -    12q1q     =   36A^B      -   36A  B    (B     +rT) 


=  -36a^b2t2  <  o     . 

We  may  conclude  that  Q'(z)  has  no  real  roots;  thus,  Q(z)  has  only  one 
real  root,  the  one  known  to  give  the  minimum  of  r(y) .   This  root  is  the 
only  local  minimum  of  r(v),  and  the  theorem  is  proved. 

The  point  at  which  the  minimum  occurs  will  be  referred  to  as 
the  pair-wise  best  point,  and  the  associated  ellipse  through  X.  and  X. 
will  be  referred  to  as  the  pair-wise  best  ellipse. 
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k.k     Three-way  Point 

If  the  functions  r(\  ,d,c2),  r(\  ,d,c2),  and  r(X  ,d,c2)  take  on 

the  same  value  at  some  point  (d,  c2)  eE,  then  there  is  an  ellipse,  a  member 

of  the  family  ^£(d,  c),  passing  through  X.,   X.,    and  \.      In  this  section 

it  will  be  shown  that  if  such  an  ellipse  exists,  it  is  unique.  An 

existence  criterion  and  parameters  of  the  ellipse  will  be  found  in 

terms  of     X.    =  x.   +iy.,    X.  =  x.  +iy.,    and  X     =  x.    +  iy    .     We  will 
111JJJ  KKK 

assume  that  0  <  x.  <  x.  <  x.  and  0  <  y.,y.,y,  • 

J.  J  K.  1   J   K 

Lemma  k . 8  If  there  exists  an  ellipse  symmetric  with  respect  to  the 

real  axis  passing  through  the  distinct  points  \.,  X.,    and  X  ,    then  it 

1    j        K. 

is  unique.   (A  general  ellipse  is  determined  by  5  points.) 

Proof  The  general  equation  for  an  ellipse  symmetric  with  respect  to 


the  real  axis  is 


2     2 

Ax  +  By  +  Cx  =  1 


Suppose  the  system 


u 


A 


Ya 


2 


A  lA\  li\ 


X. 


*k 


v     W 


has  more  than  one  solution,  say 


A. 


M 


and 


I    W 


IA  K\    M 


W  \w  VI 


then  the  vector 


IA  IA 


VI 


-1/ 


I 

7 


/ 


is  a   solution   for  any  w.      Suppose  p  /  0;    then,    for  u>  =   -B  /p  we  have 


X. 

yi 

X. 

i 

2 
x . 

2 

X  . 

*k 

2 

*k 

WA\  h\ 


w  VI 


which  implies  that 


Az      +  Cz   -   1  =   0 


has  three  distinct  roots 
implies  that 


We  can  conclude  that  ,3  =  0.   This  in  turn 


/ 


2 
x. 

l 

o 


J 


2 

y.    x. 
9±  i 

2 


M 


^k    *k 


/  W 


l°\ 

0 

W 


or  that 


az  +  yz  -   0 


8^ 

has  three  distinct  roots.  We  can  conclude  that  a  =  3  =  7  =  0,  which 

proves  the  lemma. 

Theorem  k.9     There  exists  an  ellipse  symmetric  with  respect  to  the 

real  axis  passing  through  X.,   X.,    and  K    if  and  only  if 

1   j       K 

(vxi)(yk-yi}  <  ^i^yH5  • 

If  the  ellipse  exists,  it  is  determined  by  the  parameters 


,    2,    2      2x  2,    2      2N  2,    2      2U 

d  _  1  (yi(xrxk)  +  yj(Vxi}  +  yk(xi-xj»  ( 
(y^(x  -^  +y?(vxi)  +yk(xi-xi})  ' 


a2  =  d' 


2        (yiXiXk(Xi-Xk)    +  TO^i1    +  ViVV*^ 


^.rV  +y^(vxi}  +4(vxP) 


op    a2ii    (yj(x3">)  +  yj(xfc"xi?  +y^vx3)} 

(xi-x.)(x.-xk)(Vxi) 


Proof  If  an  ellipse  is  to  pass  through  the  three  points  X.,    X.,    and  X,  , 

1     J  K. 

then  it  must  be  a  member  of  the  family  of  ellipses  passing  through  X. 

and  X,   .   Consider  the  arc  of  the  ellipse  connecting  X.    and  X,    as  we 
k  ik 

range  over  all  possible  ellipses  through  X.    and  X      (see  Figure  ^.7) • 
This  arc  sweeps  out  the  region  bounded  by  x  =  x.,  x  =  x.,  and  above  the 

±  J 

limiting  ellipse.   For  y.  ^  y,,  the  limiting  ellipse  is  a  parabola  with 

X     K 


the  equation 


x  =  Py2  +  Q , 


-  X 


PSE 


where 


and 


Figure  k.l 


Vxi 

2      2 
yk~yi 


yk-yi/ 


For  y.    =  y     it  is   the  line 
X  K. 


y  =  y±  =  yk 


The  point  X.   is  on  one  of  these  ellipses  if  and  only  if  it  is  in  this 
region,  that  is,  if  and  only  if 


x.  <  Py.  +  Q 


for   yk  -  y.  >  0  , 


x.  >  Fy  +  Q 
J     J 


for   yk  -  y±  <  0  , 


v .  >  y. 


for   yk  -  yk  =  0  . 


(See  Figure  k .8) 
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Equivalently,    we  have 


(VXi)        2     2 
(x  -x   )    <  (y.-y,)         for        y      -   y     >  0   , 

J    x       (y^-y?)      J  k       x 


(V*^        2     2 

(xj-xi)  >  7-2-2:  (yj-yi>    for   yk  -  yi  <  ° > 

(yk-yi) 


0  <  (y^)  for     yk  -  y±  =  0  , 


which  can  be  written 


(xrxi)(yk~y?)  <  (vxi)(yryi} 


for  all   cases. 

To  find  the  ellipse  parameters  d,  c2,  and  a2  for  an  ellipse 

through  \.,    \.,    and  \,    we  require  that  they  satisfy  the  three  constraint 
1   3  x 

equations: 

(a        ^2      2 
(d-x. )     y 

}  a2    +  a2  -  c2  =   ' 

(a         \2  2 

(d-x  )     y. 

pi    a +  a —  -  1 

y      a2     a2  -  c2     ' 

2      2 
(d-Xk)      yk 
^j      a2    +  a2  - c2  =    " 

Assume  that  \. ,    \.,    and  \_  satisfy  the  existence  criterion.   Then,  \.   must 
lie  in  the  region  described  above  (see  Figure  ^.8).  We  can  establish 
three  cases: 


Case      .  Cf  yk-y.  ,    then  y  .-.v. 

Case   II.         It'  Yk~yi   <  !    n  yi"yK  ' 

ie  in.     If  yk-yi  ■  0,    then  y,-yk 

y.-y.    >  0. 


*i 


•Xi 


CASE    I 


CASE  n 


CASE  HI 


Figure  4.8 

Case   I.      We  know  from  Section  k.J   that   if  y.    -y.    4  0,    we   can  use 

Jk   i 

constraint  equations  1)  and  5)  to  get  an  expression  for  a2  in  terms 
of  d.   We  have 


a2  =  (d  -  ( 


\ 


+x, 


Vxi  yk+yis , , .      A+Xi     Vxi  yk"yi 


2   y,  -y. 

^k  Ji 


■))(d  - 


yi  +y- 

Jk  '  i 


)) 


(d-P13)(d-Q13) 


Likewise,  since  y.  -y.  /  0,  we  can  use  constraint  equations  l)  and  2) 
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to  get 

x.+x.       x.-x.   y.+y.  x.+x.       x.-x.   y.-y. 

a2=  (d-  (^-V^(d-   ^"^y^ 

=    (d-P12)(d  -Q^)    . 

Setting  the  two  right-hand   sides  equal  and  solving   for  d,    we  get 
(d  -P13)(d  -Q15)    =    (d  -P12)(d  -Q^)    , 


and 


^3^3  "  P12Q12 


Plugging  in  the  proper  expressions,  we  get 


2,    2      2,  2/2      2v  2,   2      2, 

d_   lyj^j-V    +yj(VXi}    +yk(VXj) 

"  y^x  -x^  +  yj(xk-x.)  +  yk(x.-x.) 


(Notice  that  a  positive  denominator  is  equivalent  to  the  existence 

criterion. ) 

If  we  write  d  =  —  —  and  plug  this  back  into  the   equation 
d   K 

a2  =    (d-P13)(d-Q15)    =  d2   -    (P13  +Q13)d+P13Q15   , 


we  set 


a2  .  a2  -  I  (p15  +Ql5)  I  +  P^ 

^(P13+QI3)j-gP13Q13K 

2K 


the  proper  expressions    ind   rearrani  Li  ms,   wo  gi 


a2  m     i    _  fojVV*^    '  ^V*!* 


y?(xfV   '  yj(W   '  yk'    l-xj5 


From  Section  U o  we  know  that  since  y  /  y.  we  can  also  use  constraint 

K  1 

equations   l)    and  3)    to  get 

2     2 

Vxi         i    yk"yi 
(d-(JL-i  +  i^— i)) 
2  2   X-x. 

c2  =   a2  L_i_ 

X,  +x. 

(yk-yi} 

a2  f  1 — — 


(2d(VX.).(^)) 


Again  using  d  =  —  — ,    we  have 


1  J 

2  K' 


(yk"yi)K 
c2  =   a2   I  1 -K 


((xk-xi)j  -(x^-x^)K) 


which  yields 


c2-a2fl       ^V**'   IJ^VV    +  yk<Vxj! 

(x.-x^lx.-x^K^-x.) 


Case  II.   In  this  case  y  -y.  ^   0  and  y .  -  y  ^'  0,  so  we  may  use  constraint 

K   l  j   k 

equations  1)  and  3),  and.  2)  and  3)  to  get  expressions 

a2  =  (d  -P13)(d  -C^,)  , 
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and 

a2  =  (d-P23)(d-Q23)  . 

Using  the  same  procedure  as  in  Case  I,  we  arrive  at  the  same  expressions 
for  d,  a2,  and  c2. 

Case  III.   In  this  case  y.  -y.  ^  0  and  y.  -y  ^   0.  Using  constraint 
equations  I)  and  2),  and  2)  and  3),  we  again  achieve  the  same  results. 
This  proves  the  theorem. 

The  point  at  which  three  functions  take  on  the  same  value  will 
be  referred  to  as  a  three-way  point,  and  the  associated  ellipse  through 
X.,    X.,    and  X     will  be  referred  to  as  a  three-way  ellipse. 

1    J         K 

k*5     The  Algorithm 

The  results  of  the  first  four  sections  of  this  chapter  have 
shown  that  the  solution  of  the  mini -max  problem, 


min     max  r(\.,d,  c2)  , 
(d,c2)eR  X.eE+         X 


is  the  minimum  point  of  a  single  function,  a  pairwise  best  point,  or  a 
three-way  point.   In  this  section  an  algorithm  will  be  developed  to  find 
the  mini -max  solution  among  the  possible  candidates. 

From  the  Alternative  Theorem  we  know  that  if  a  point  is  the 
mini -max  solution,  the  function  or  functions  it  is  associated  with  must 
equal  the  maximum  over  all  functions  at  that  point .  Equivalently,  the 
associated  ellipse  must  contain  the  positive  hull,  H+,  in  the  closure 
of  its  interior.   For  this  reason  the  minimum  point  of  a  single  function 
cannot  be  the  mini -max  solution. 


l.omina   h  .  1  '  tie    pOfi  luil  1  , 

ilue,    then  the  mini -max   solution  cannot  be  the  minimum 

le    function. 
Krom  the  Alternative  Theorem  we  know  that  the  mini -max  so 
occur  at   a   local  minimum  of  a   ^:in;rle   function,    say  r(\.,d, ci     , 

J 


on]y  if 


r(\^,d,  c2)   :      maxi    r(\. ,  d,  c2) 

L.  i 

I 


'    •>      '  +      v   i- 

X.ell 


at  that  point.      That  is   to   say,    the  associated  ellipse  through  X. 
contains   the  positive  hull  in  the   closure  of  its   interior.      From 

Theorem  k  .3  we  know  that  the  only  local  minimum  of  r(\.,d,  c2)   occurs 

J 

at  (d, c2)  =  (x.,-y.).   The  associated  ellipse  is  the  degenerate  ellipse 

J      J 

with  foci  at  X.    and  X..      If  H  contains  more  than  one  eigenvalue,  say 

J  J 

K    ^  X.,    then  K    must  lie  outside  the  degenerate  ellipse.   Thus,  we  have 

r(W-yj>  >  r(X.,x.,-y^  , 

which  proves  the  lemma. 

We  now  know  that,  the  mini -max  solution  is  either  a  pair-wise 

best  point  or  a  three-way  point. 

Tlieorem  l4-.ll  A  pair-wise  best  point  whose  associated  ellipse  contains 

the  positive  hull,  H+,  in  the  closure  of  its  interior  is  the  mini-max 

solution.   If  no  such  pair-wise  best  point  exists,  then  the  mini-max 

solution  can  be  found  among  the  three-way  points  whose  associated  ellipses 

contain  the  positive  hull. 

Proof  Suppose  the  pair-wise  best  point  associated  with  X.   and  X   ,  call 

it  (d.,  ,c2.,  ),  is  such  that  the  associated  ellipse  passing  through  X. 
'  jk'   ok  J 

and  \,  contains  H+.   Then,  we  have 

k  ' 
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r(WcV=  ™*+r(WcV 


min    max  r(\.,d,c2) 
(d,  c2)eR  A..eH+    X 


min    max{r(\.,d,c2),r(\,,d,c2)] 
(d,c2)eR        3  K 


=  r(VVc2^  ? 


which  proves  the  first  result.   The  last  equality  is  a  result  of 
Lemma  J+.10  and  the  Alternative  Theorem.   The  mini -max  problem  over  two 
functions  must  have  a  pair-wise  best  point  for  its  solution. 

Suppose  that  no  pair-wise  best  point  qualifies.   From  the 
Alternative  Theorem  and  Lemma  4.10,  we  know  that  the  mini -max  solution 
is  a  three-way  point  whose  associated  ellipse  contains  the  positive  hull. 
This  proves  the  theorem. 

To  find  the  mini-max  solution,  one  must  systematically  take 
each  pair  of  eigenvalues  in  the  positive  hull  and  find  the  pair-wise  best 
point.   If  the  associated  ellipse  contains  the  positive  hull,  then  the 
point  is  the  mini-max  solution.   If  no  pair-wise  best  point  qualifies, 
one  must  take  each  combination  of  three  eigenvalues  and  look  for  a  three- 
way  point.   If  it  exists  and  its  associated  ellipse  contains  the  positive 
hull,  the  point  is  a  candidate.   The  mini-max  solution  is  the  three-way 
candidate  that  yields  the  smallest  convergence  factor,  r(\, d, c2) . 

Notice  that  those  eigenvalues  in  H4"  which  are  involved  in 
some  combination  of  three  eigenvalues  that  produced  a  three-way  candidate 


are  exactly  the  key  elements  described  in  Section  ..  .   [n  bl     '-3e 

be  mini -max  solution,  the  key  «      t,s  are  determined. 
Thus,  we  may  discard  those  eigenvalues  in  H"f  which  are  not  key  elements 
Reducing  the  number  of  eigenvalues  in  If'  redu     tie  number- 
combinations  of  three  that  must  be  tried  in  subsequent  searches  for 
a  mini -max  solution. 
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5.   ADAPTIVE  PROCEDURE 

The  algorithm  for  finding  optimal  iteration  parameters 
developed  in  Chapter  2  and  Chapter  3  depends  upon  knowledge  of  the 
spectrum  of  the  matrix  A.   In  practice,  however,  little  is  known  about 
the  spectrum  of  A.   In  this  chapter  a  procedure  will  be  developed  for 
estimating  the  hull  of  the  spectrum  from  information  acquired  during 
the  iteration,  allowing  dynamic  improvement  of  the  iteration  parameters. 
Section  5.1  will  show  how  eigenvalue  estimates  can  be  extracted  from  the 
sequence  of  residuals.   Section  5-2  will  show  how  these  estimates  can 
be  used  to  build  an  approximate  hull  of  the  eigenvalues  of  A. 

5-1  Modified  Power  Method 

Suppose  that  the  iteration  parameters  d  and  c2  have  been  chosen 
from  some  prior  knowledge  of  the  matrix  A.  After  n  steps  of  the  Stiefel 
Iteration  based  upon  d  and  e2,  the  error  can  be  expressed  as 


e  =  P  (A)e_  . 
n    n    0 


Multiplying  this  expression  by  A,  we  can  write  the  residual  as 


r  -  P  (A)rn  . 
n    n    0 


Suppose  that  r  can  be  written  as  a  linear  combination  of  the 
eigenvectors  of  A;  then,  we  have 


r_  =  Z  (a. v.  +<xv. )  , 
0   .  _  v  1  1   1  i'  ' 

i=l 


wlu  ■  Ls   tin-  i'i  ;envector  of  A  associated   with   the  eigei  \ 

and    ||v.  ||  =   1.      Since  A  is  a  real  valued  matrix  and  r      I 
valued  vector,    the  eigenvalues   and  eigenvectors  appear-    ir 
conjugate  pairs.      After  n  steps  of*  the   iteration,   we  have 


r     =   P   (A)r,  =     l\      («.  P   (\.)v.   +a.P   (T..)v.) 
n         nv  .    .    v  i  n     i'   i       1  n^   i'    iy 

1=] 


From  Section  3.1  we  know  that 


Pn(\)  =  (M(^))n 


for  large  n.   Let 

m.  =  M(\. ) 
l      l 

be  the  eigenvalue  of  the  operator  M(A)    corresponding   to   the  eigenvalue 
\.    of  A.      For    Large   n  we   can  write 


•     .,  /       n         — u — 

r     =     £  (a.m. v.   +a.m.v.  ) 
n       .    ..        i  i  i       ill7 

i=l 


Suppose  m,   and  m     are  the  dominant  and   subdominant  eigenvalue: 

Q.  S 

of  M(A) ;    i.e.,    suppose 

|m,  I   >    |m    |   >    jm.  |      for     i  J  d,  s    . 

1    d     —    '    s  '  i  '  '       ' 

Then,    for  large  n  the   residual  can  be  written  as 

n,  v     _n, v        n,  v     _n, N 

r     =  m   (a,v, )  +m,  (a,v, )  +m  (a  v   )  +m   (a  v   )  +€     , 
n         add         dv   d  d  sv    s   s  sv    s   s  n 

where  c  is  small  relative  to  the  other  terms.   The  residual  is  nearly 
n  J 

a  linear  combination  of  the  eigenvectors  associated  with  the  dominant 
convergence  factors.   Likewise,  we  have 
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n+1,  v     —n+1, v        n+1,  N     —n+1, >, 

r    •    =  m,      (a,  v..)  +m,      (a,v,)  +m       (a  v  )  +m       (a  v  )  +e     . 
n+1         d         dd'        d      vddy        s  s   s'        s  ssy        n+1 


n+2,  s      -n+2, \        n+2,  N     -n+2, v 

r         =  m       (a,vj  +m,      (a,v,)  +m       (a  v   )  +m       (a  v  )  +e     „  , 

n+2         d         dd         d         dd         s  s   s  s  s   s         n+2  ' 


n+3  /  n     _n+3  / \       n+3  /  \     _n+3  / \ 

r         =  m ,      (a,v,)  +m,      (a,v,)  +m       (a  v  )  +m       (a  v  )  +e     ,   . 

n+3         d         d  d'        d         dd         s         ss  s         ss         n+3 

n+4,  v     _n+4, s       n+4,  N     _n+4, N 

r    ,,,  =  m,      (<x,v,)  +m..      (a,v,)  +m       (a  v   )  +m       (a  v  )  +e     .     . 
n+4         d         dd'        d         dd'        s      vss/        s      v    s   s'        n+4 


If  we  let 


h  3  2 

Q(z)    =    z     +qjzJ  +(^2     +  q^z  +  q^ 


(z-m,)  (z-m, )  (z-m  )(z-m  )    . 
v       d'  v       d  s  b     ' 


then,    ignoring  the   e  terms,    we  have 

III"     i    +Q.-.r     -z  +  cur     ^  +  q*r     i   +Q.i  r    II   =   0    . 
11  n+4     Hl  n+3       *?  n+2     ^3  n+1     H4  n" 

Choosing  (q  ,  q~,  q,,  q,  )  to  make  the  I  -norm  of  the  above  expression  as 
small  as  possible  is  equivalent  to  finding  the  least  squares  solution 
of  the   system 


n+3       n+2       n+1       n 


The   roots  of  the  polynomial 


4  3  2 

Q(z)    =    z     +q1z/+q2z     +q^z+q^   , 


are   bh<  ition  o 

abovo,    are  approximations   bo  the  dominanl 

In  Light  of  the  asymptotic  equivalence  of   the   two  operate 

Pn(A)    ~    (M(A))n   , 

this  procedure  corresponds  to  the  power  method  for  simultaneous  estimation 
of  the  eigenvalues  of  M(A)  (Wilkinson  [2^]).   The  number  of  estimates 
produced  is  determined  by  the  degree  of  the  polynomial  Q(z)  and  the 
number  of  previous  residual  vectors  stored.   The  method  can  be  adjusted 
to  yield  any  number  of  eigenvalue  estimates. 

The  accuracy  of  this  approximation  depends  upon  the  separation 
of  the  eigenvalues  of  M(A)  (Wilkinson  [2h])  .      For  the  purposes  here,  it 
is  only  important  that  the  error  lies  in  the  direction  of  the  other 
eigenvalues  of  M(A) .   If  A  is  symmetric,  this  will  hold. 
Theorem  5 ■ 1  If  A  is  symmetric,  the  eigenvalue  approximations  produced 
from  the  procedure  described  above  will  lie  in  the  convex  hull  of  the 
eigenvalues  of  M(A) . 

Proof  If  A  is  symmetric,  then  it  has  a  complete  orthonormal  set  of 
eigenvectors.   If 

Q(z)  =  z  -i-q  z  4-q^z  +  q  z  +  q^ 


(Z-Tl1)'(z-Tl1)(z-T]2)(z-Tf2)  , 


then  for  large  n  we  have 
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'     ,    +  qn  r     ,  +  q_r      _  +  q_,r     _   +  q,  r 
n+k     Hl  n+3     ^  n+2     H3  n+1     ^k  n1 


i=l 


Z      |a.|2|mi|2n|Q(m.)|2 
i=l 


If    |  Q,(m. )  |    can  be  made   smaller  for  each  m. ,    then  the  entire   sum  will 
be   smaller.     We  have 

IqOil)  I   =    |mi-T]1|  |mi-Tf1|  ^-Tigl  |mi-n2l    • 

Suppose  t\     lies  outside  the  convex  hull  of  the  eigenvalues  of  M(A) 
(see  Figure  5.1).      Both    |m.  -tj    |    and    Im.-rj"    |    can  be  made   smaller  for 


t 

{ 

< 

• 

• 

• 

• 

i       \ 

Figure   5*1 


;.  run  |  e  convex  hul 

l)  .   By  tin  'urn  can  b<         by  mo\ 

toward  the  convex  hull.  We  may  cone1  .    hat  the  Leasl 

ution  will  yield  eigenvalue  approximations  that  lie  in      he 
convex  hull,  which  proves  the  theorem. 

No  such  result  is  known  by  the  author  for  a  general  matrix, 
^ed,  the  loss  of  orthogonality  and  the  possibility  of  nonlinear 
elementary  divisors  complicate  the  analysis.   Experimental  results 
on  large  sparse  matrices,  however,  have  exhibited  the  tendency  of  the 
approximate  eigenvalue  to  lie  inside  the  convex  hull  of  the  true 
eigenvalues . 

The  relationship  between  eigenvalues  of  A  and  eigenvalues  of 
M(A)  is  given  by  the  function 

(cosh~"(— ^— )  -cosh"  (— )) 
M(X)  =  e        C  C 

j_ 

(d-X)  +((d->Q2  -c2)2 
1 

d  +(d2-c2)2 

This  function  maps  the  members  of  the  family  *£(d,  c)  onto  circles 
centered  at  the  origin  (Section  2.1).   The  radius  of  the  circle  is 

!m(\)  j  =  r(\,d,  c2)  , 

the  convergence  factor  associated  with  the  ellipse.   The  region  of 
convergence  in  the  \-plane,  shown  in  Figure  5-2,  is  mapped  injectively 
onto  the  region  in  the  m-plane  shown  in  Figure  5.2.   Since  the  branch 
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X-PLANE 


m- PLANE 


Figure   5.2 


of  cosh"  "  with  positive  real  part  is  used,  we  have 


|M(\) 


cosh"  (— :— )  -cosh"  (-) ) 

\  c  VC 


>  e 


■cosh  (-) 


<§)♦  ((f)2  -ifi 


c2 


d  +(d-c2) 


Let 


g  =  d+(d2-c2)2   ; 


thou,  we  can  write 


with  the  restriction 


^  =  d  -  -  (^  (  -)  , 


-  g 


If  poor  separation  of  the  eigenvalues  of  M(A)  causes  the  modified  power 
method  to  yield  an  eigenvalue  approximation  m  such  that 


Iml  <  |C2 


g 

then  it  must  be  discarded. 

This  map  is  conformal,  but  it  is  not  linear.   If  m  lies 
in  the  convex  hull  of  [m.  ] ,  X   does  not  necessarily  lie  in 
the  convex  hull  of  [X. } .     An  underestimation  of  the  imaginary  part  of 
m.  may  cause  an  over estimation  of  the  real  part  of  d-\. .   This  warping  is 
slight,  however,  and  shows  little  effect  in  practice. 


5.2  Procedure  and  Example 

A  four  step  procedure  for  dynamic  improvement  of  iteration 
parameters  is  outlined  below. 
1.   Choose  initial  d  and  c2. 

The  initial  choice  of  d  and  c2  must  be  based  upon  prior 
knowledge  of  the  matrix  A.   In  Section  1.2,  crude  bounds  were  established 
which  gave  rectangular  regions  known  to  contain  the  spectrum  of  A.   These 
regions  could  be  used  to  choose  initial  d  and  c2.   If  the  system  were 
scaled  so  that  each  diagonal  term  of  A  had  the  value  1,  then  d  =  1, 
c2  =  0,  representing  circles  centered  at  1,  could  be  used. 
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2.   Iterate. 

Perform  a  predetermined  number  of  steps  of  the  Stiefel 
iteration  based  upon  d  and  c2,  and  store  the  last  five  residuals. 
The  number  of  steps  should  be  large  enough  to  insure  separation  for 
the  modified  power  method. 
3-   Get  eigenvalue  estimates. 

Using  the  modified  power  method  described  in  Section  5«1> 
obtain  eigenvalue  estimates.  Add  the  new  eigenvalue  estimates  to  any 
previous  eigenvalue  estimates  and  form  the  positive  hull. 
k.      Choose  new  d  and  c2. 

Using  the  algorithm  outlined  in  Section  h.^,    find  the  optimal 
d  and  c2  for  this  set  of  eigenvalue  estimates. 

Repeat  2,  3,  and  k.      These  steps  will  be  referred  to  as  a 
cycle. 

To  illustrate  how  this  procedure  might  work,  suppose  the  hull 
of  the  eigenvalues  of  A  is  as  shown  in  Figure  5-3 •   Suppose  d  and  c2 
have  been  chosen  as  indicated  in  Figure  ^.h.      Several  members  of  ^(d,  c) 
are  also  shown  in  Figure  ^.k.      It  is  clear  that  X     and  X     are  the  dominant 
and  subdominant  eigenvalues  for  this  choice  of  d  and  c2. 

After  a  sufficient  number  of  steps  of  the  Stiefel  iteration, 
the  modified  power  method  will  give  eigenvalue  estimates  p  and  p  shown 
in  Figure  5«5«   If  we  let  d+c  and  d-c  also  be  eigenvalue  estimates,  then 
the  approximate  hull  is  as  shown  in  Figure  5-5.  A  choice  of  d  and  c2 
based  upon  this  approximate  hull  yields  the  d  and  c  shown  in  Figure  5-6. 
Some  members  of  the  new  family  *%-(&,  c)  are  shown  in  Figure  5-6,  and  it 


The  eigenvalue  estimates  were  chosen  for  illustrative  purposes.   The  values 
of  d  and  c2  were  computed  from  these  values,  and  the  figures  were  plotted 
by  the  IBM  CALCOMP  Plotter  at  the  University  of  Illinois. 
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Figure   5 '5 


Figure  5.6 


is  clear  that  Xr   and  \f  are  the  dominant  and       inant  eigenvalues 
for  this  choice  of  d  and  c2. 

After  a  sufficient  number  of  steps  of  the  Gtiefel  iteration, 
the  modified  power  method  will  give  eigenvalue  approximations  p  an<i 
as  shown  in  Figure  5 •! •      The  new  approximate  hull,  shown  in  Figure  5-1, 
describes  the  true  hull  fairly  well.   The  optimal  choice  of  d  and  c2 
based  upon  this  approximate  hull  yields  the  d  and  c  shown  in  Figure  5.8. 
The  best  ellipse  enclosing  the  approximate  hull  is  labeled  "computed"  in 
Figure  5-8,.  and  the  best  ellipse  enclosing  the  true  hull  is  labeled 
"optimal".   They  differ  only  slightly. 

If  further  eigenvalue  estimates  fall  inside  the  approximate 
hull,  they  will  provide  no  new  information  and  will  be  ignored.   If 
further  eigenvalue  estimates  lie  outside  the  approximate  hull  but  inside 
the  true  hull,  they  will  help  to  better  describe  the  true  hull  and  will 
lead  to  better  values  of  d  and  c2. 

Notice  that  in  this  example  the  true  hull  was  inside  the  region 
of  convergence,  the  ellipse  passing  through  the  origin,  for  each  choice  of 
d  and  c2  (see  Figure  ^.K).      This  may  not  always  be  the  case.   If  an 
eigenvalue  X.    is  outside  of  the  region  of  convergence,  then  the 
corresponding  eigenvalue  m.  of  M(A)  will  have  modulous  greater  than 
unity.   The  error  in  the  direction  of  the  eigenvector  associated  with  X. 
will  increase  at  each  step.  However,  m.  will  dominate  so  strongly  that 
at  the  end  of  the  cycle,  the  modified  power  method  will  yield  an  excellent 
approximation  of  X.  .   The  next  choice  of  d  and  c2  will  involve  this 
estimate.   If  a  provision  is  made  to  return  to  the  iterate  at  the  beginning 
of  the  cycle,  a  poor  choice  of  d  and  c2  will  not  move  the  iteration  toward 
convergence,  but  it  will  give  rise  to  new  iteration  parameters  that  will 
produce  convergence . 
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Figure  5 • 7 


OPTIMAL 


Figure   5-8 


6.   IMPLEMENTATION  AND  RESULTS 

The  implementation  of  the  Tchebychef  algorithm  will  be 
discussed  in  Section  6.1,  followed  by  a  discussion  of  competing  methods 
in  Section  6.2,  and  experimental  results  in  Section  6.3.  A  listing  of 
the  subroutines  coded  in  FORTRAN  will  appear  in  the  appendix. 

6 . 1  Implementation 

The  Stiefel  iteration  is  performed  by  the  subroutine  TCHEB 
(Appendix  A).  A  basic  flowchart  of  TCHEB  appears  in  Figure  6.1.   To 
solve  the  system  Ax  =  b,  the  user  must  supply  the  matrix  A  of  dimension 
N,  a  storage  scheme,  and  a  subroutine  to  perform  matrix  vector 
multiplication  (NSYMAX  in  the  listing) .   The  user  must  also  supply  an 
initial  guess  at  the  solution  x  (x(N))  and  the  target  vector  b  (B(N)). 
Storage  space  must  be  provided  for: 

R(N)  the  residual, 

DX(N)  the  change  in  X  at  each  step, 

XLAST(N)  X  at  the  start  of  the  cycle, 

RLAST(N)*  R  at  the  start  of  the  cycle, 

S(N,U)       the  four  previous  residuals  at  the  end 
of  the  cycle. 

The  dimension  parameter  N  must  be  passed  in  a  common  statement 
as  well  as  the  parameter  ND,  the  number  of  diagonals  required  to  store  A, 


This  vector  may  be  replaced  by  a  matrix  vector  multiplication  in  the 
subroutine  LASTX  (Appendix  A) . 
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DX 
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ADAPTIVE     PROCEDURE 
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THEN     RESTART    STIEFEL 
ELSE     RESUME    STIEFEL 


Figure  6 . 1 


The  initial  choice  of  d  and  c:  must  also  be  par..     i  a 
common  statement.   This  choice  can  be  based  on  rough  bo  ■ 
those  established  in  Section  1.2.  A  poor  choice  of  d  and  c2  will 
cause  the  first  cycle  to  be  spent  in  obtaining  eigenvalue  estima1 
(Section  5 .2).   Care  should  be  taken  to  avoid  a  choice  that  causes 
overflow  before  the  end  of  the  cycle,  at  which  time  the  program  will 
adjust  itself. 

The  positive  hull  of  the  eigenvalue  estimates  is  stored  in  the 
doubly  linked  list  CH(25,2)  with  link  array  ICH(25,2).   The  foci  d+c 
and  d-c  associated  with  the  initial  choice  of  d  and  c2  are  used  as  initial 
eigenvalue  estimates.   For  that  reason  d  and  c2  should  be  chosen  so  that 
d+c  and  d-c  lie  inside  the  convex  hull  of  the  eigenvalues  of  A.   If  the 
user  wishes  to  supply  a  priori  eigenvalues,  he  must  initialize  CH  and 
ICH  in  the  subroutine  IWIT  (Appendix  A) . 

The  number  of  iterative  steps  per  cycle  must  be  set  by  the 
input  parameter  ICYCLE.   In  practice,  little  difference  was  found  for 
values  of  ICYCLE  from  15  to  50.   The  value  20  was  used  in  most  of  the 
results  that  follow. 

The  other  input  parameters  are  IMAX,  the  maximum  number  of 
iterations  allowed,  and  EKBND,  the  acceptable  bound  on  the  i   -norm 
of  the  residual  vector.   Since  the  hull  of  the  eigenvalues  of  A  is  found 
by  the  adaptive  procedure,  a  bound  on  the  modulus  of  the  smallest 
eigenvalue  is  available.   This  may  be  incorporated  in  the  error  bound  by 
replacing  EKBND  by  EKBND* (D-DSQRT(A2))  in  the  test  to  halt  in  subroutine 
TCHEB. 

The  adaptive  procedure,  called  from  subroutine  TCHEB,  is 
broken  into  four  parts  (see  subroutine  ADAPT,  Appendix  A).   The  first 
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part,  subroutine  LSTSQ,,  involves'  finding  the  least  squares  solution 
to  the  system 

where  the  columns  of  S  are  the  previous  residuals  (Section  5.1). 
Subroutine  LSTSQ  converts  this  system  to  the  equivalent  system 


T  T 

S  S,   i  *  Q,   n  =  -  S  R   .  . 

4X4     4X1         4X1 


Since  the  residuals  may  become  small,  this  system  is  multiplied  by  a 

T  T 

constant  so  that  S  S(l, l)  =  1.0.   The  matrix  S  S  is  symmetric,  so  it  can 

be  computed  by  taking  10  inner  products.   Likewise,  4  inner  products  are 

T 
required  to  compute  -S  R.   The  normalized  system  is,  in  general,  less 

stable  than  the  large  system,  but  results  have  shown  that  this  instability 

does  not  warrant  the  extra  work  involved  in  solving  the  large  system. 

The  normalized  system  may  be  solved  in  many  ways.   If  one 

T 
eigenvalue  of  M(A)  dominates  strongly,  however,  the  matrix  S  S  may  be 

singular  (Section  5.1),  in  which  case  a  least  squares  solution  is  desired. 

For  that  reason  a  Bidiagonalization  method  with  reorthogonalization  was 

used  to  solve  the  4  x4  system  (Golub  and  Kahan  [11]).   The  subroutines 

for  the  Bidiagonalization  method  appear  in  Appendix  B. 

Subroutine  EVS  finds  the  roots  of  the  polynomial  associated 

with  the  least  squares  solution  Q  and  converts  them  into  eigenvalue 

estimates  (Section  5*1) •  Any  root  finding  routine  may  be  used  to  find 

the  roots  of  the  fourth  degree  polynomial.  A  library  routine,  RSSR, 

from  the  University  of  Illinois  F0RTU0I  Library  was  used  by  the  author. 

A  listing  of  RSSR  appears  in  Appendix  C. 


The  subroutine  HULL  adds  the  new  eigenvalue  estimai 
the  previous  eigenvalue  estimates  and  forms  the  positive  hull  (Sectio.-. 
[f  none  of  the  new  eigenvalue  estimates  are  members  of  the  po 
hull,  control  is  returned  to  the  subroutine  TCIIEB,  and  the  iteration 
is  resumed. 

The  subroutine  ELLIP  finds  the  best  ellipse  enclosing  the 
positive  hull  (Section  ^.5)«   If  the  best  ellipse  is  a  three-way  ellipse, 
this  subroutine  also  determines  the  key  elements  of  the  positive  hull  and 
discards  the  others  (Section  3-2  and  Section  k.k).      Finding  a  pair-wise 
best  ellipse  involves  finding  the  root  of  a  fifth  degree  polynomial  in  a 
certain  interval  (Section  h.J>)  .   This  was  done  by  the  routine  ZEROIN, 
which  appears  in  Shampine  and  Allen  [l8].  A  version  adapted  for  the 
purposes  here  appears  in  Appendix  C. 

6.2  Competing  Algorithms 

Two  competing  iterative  methods  for  solving  the  nonsymmetric 
linear  system  Ax  =  b  are  the  Bidiagonalization  method  (Golub  and  Kahan 

[11],  Paige  [17])  and  the  method  of  Conjugate  Gradients  applied  to  the 

T      T 
equivalent  system,  A  Ax  =  A  b  (Hestenes  and  Stiefel  [12]).   The  Tchebychef 

algorithm  requires  more  storage  than  the  other  two  in  that  four  previous 

residuals  must  be  retained  for  the  adaptive  procedure.   The  other  methods 

require  more  work  per  iterative  step.   Table  6.1  shows  a  comparison  of 

the  work  per  iterative  step  and  storage  requirements  of  the  three  methods. 
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Table 

6.1 

mul/step 

add/step 

STORAGE 

TCHEB 

(ND+2.7)N* 

(nd+3-7)n 

(ND+9)N 

BIDIAG 

(2ND+7)N 

(2ND+5)N 

(nd+5)n 

CG  on  ATA 

(2ND+6)N 

(2ND+6)N 

(ND+5)N 

Hero  N  is  the  dimension  of  the  system,  and  ND  is  the  number 
of  columns  of  length  N  required  to  store  the  matrix  A.   If  A  is  a 
5-point  difference  operator,  then  ND  =  5*   Finite  element  matrices 
usually  require  more  storage.   For  instance,  when  cubic  shape  functions 
are  used  in  two  dimensions,  we  may  have  ND  =  33  and  in  three  dimensions 
we  may  have  ND  =  135  (Zienkiewicz  [27]) •  As  the  storage  required  for 
the  matrix  A  becomes  greater,  the  extra  storage  required  by  the  Tchebychef 
algorithm  becomes  less  of  a  factor,  but  the  work  required  by  the  other 
methods  remains  twice  as  much  as  that  required  by  the  Tchebychef  algorithm. 

Both  the  Bidiagonalization  method  and  the  method  of  Conjugate 

T      T 
Gradients  on  the  equivalent  system,  A  Ax  =  A  b,  are  sensitive  to  the 

T 
condition  of  A  A,  whereas  the  Tchebychef  method  is  sensitive  to  the 

condition  of  A.   Because  of  this  advantage  and  the  advantage  of  less 

work  per  iterative  step,  the  Tchebychef  method  was  considerably  faster 

than  the  other  methods  on  a  series  of  test  problems.   In  all  of  the  tests 

that  were  run,  the  Bidiagonalization  method  was  slightly  better  than  the 

method  of  Conjugate  Gradients.  For  that  reason  results  are  shown  for 

the  Bidiagonalization  method  only. 


The  adaptive  procedure  in  the  Tchebychef  algorithm  requires  l^N 
multiplications  and  additions  (Section  6.1).   If  an  adaptive  procedure 
is  carried  out  every  20  steps  (ICYCLE  =  20)  then  the  average  overhead 
is  .7N  per  step. 


Results 

The  convergence  properties  of  the  Tehebychef  algoritlu     nd 
only  upon  the  spectrum  of  the  matrix  A  and  not  upon  any  special  zero 
structure.   It  is  desirable,  however,  to  test  the  algorithm  on  easily 
constructable  systems  whose  spectrums  have  known  properties.  Consider 
the  second  order  linear  differential  operator  on  a  a  rectangular  domain: 

-  |-  (A  (x,y)  A)  .  ±  (A  (x,y)  A)  +  B_(x,y)  A  +  Bp(x,y)  |-  +  Q(x,y)  , 
dx   J-     dx    dy  d  dy     L  dx    d  dy 

where  the  functions  A  ,    Ap,    B, ,    Bp,    and  Q  are  continuous   and  dif ferentiable 
This  operator  can  be  broken  into  two  parts: 

-  |-   (A    (x,y)   ±)    -  A.   (A   (x,y)   f )    +  Q(x,y)    , 

dx       L  dx         dy       d  dy 


and 


yx,y)  I  +  B2(x,y)  A. 


With  suitable  boundary  conditions,  the  first  part  is  a  self-adjoint 
operator  while  the  second  is  not  (Dunford  and  Schwartz  [^]).  Accordingly, 
if  we  approximate  this  operator  by  a  finite  difference  operator  on  a 
regular  grid,  the  first  part  gives  rise  to  a  symmetric  matrix,  and  the 
second  part  gives  rise  to  a  nonsymmetric  matrix.   Let  h  be  the  grid  size 
and  (x.,y.)  be  the  grid  points  with  the  standard  left-to-right,  down-to-up 
ordering  of  node  points.  Using  central  differences  and  Dirichlet  boundary 
conditions,  let  A  be  the  5-point  difference  matrix  associated  with  the 
differential  operator.   The  matrix  A  can  be  written  in  two  parts: 
A  =  M+N,  corresponding  to  the  two  parts  of  the  differential  operator. 
The  matrix  M  is  symmetric  and,  if  k  is  the  number  of  grid  points  per 


column,    can  "be  written  in  block  form  as 
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M 


\ 


1 

M2, 


'k-1 


'k-1 


\ 


The  M. 's  are  tridiagonal  and  the  C. 's  are  diagonal.   If  l   is  the  number 
of  gridpoints  per  row,  we  have 


and 


M.  = 

1 


li 


2i     2i 


'2i    "3i 


"  b 


£-li 


'  Vii  an 


'2i 


n 


where 


ai,j  -7?  lAi(xi  "l'3^  "  \(xi  4  ?'yj} 


+A2(x.,y.    -§)    +  A2(x.,y.    +  §) }    +  Q(x.,y.)    , 


1a/  h  x 


1  A        /  h\ 

Cij=^2A2(Vyj    +  2> 


The  matrix  N  is  nonsymmetric  and  can  be  written  in  block 


form  as 


N  = 


N     F 

1  1 

F     N  F 

2  2  2 

-F  N 

3  3 


k-l 


■  F.    N. 
k    k 


Again,  the  N. ' s  are  tridiagonal  and  the  F. 's  are  diagonal.  We  have 


N. 


r 

ali 

1  d£1 

0 

d21 

3i 

0 

"£-li 


*-d0.    0 


n6 


and 


F.  = 

l 


"2i 


where 


d.  .  =  -prr-   B-  (x.,y.)  , 


fiJ  "  2h  %^'j) 


a)   Nonsingular  5-point  Difference  Matrices 

If  A  (x,y),  A  (x,y)  >  0  and  Q(x,y)  >  0,  then  M  is  a  positive 
definite  matrix.   If  B-  (x,y)  and  Bp(x,y)  are  constant  functions,  then 
N  is  antisymmetric  (Varga  [22]).   From  Section  1.2  we  know  that  the 
eigenvalues  of  A  =  M+N  lie  in  the  right  half  plane.   In  particular,  let 


A1(x,y)  =  Ag(x,y)  =  1  , 


Q(x,y)  =  0  , 

B1(x,y)  =  B2(x,y)  =  P  , 

on  the  square  region 

0  <  x  <  41,   0  <  y  <  kl   , 

with  grid  length  h  =  1.   Then,  A  is  a  1600  xl600  matrix  and 


L17 


M,    = 


-1 


'    .    0 


'-1 
■  1    '   k 


c. 

1 


-1 


-1 


N. 

l 


2 


F.    = 
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The  eigenvalues  of  the  positive  definite  matrix  M  lie  in 
the  open  interval  (0,8),  while  the  eigenvalues  of  the  antisymmetric 
matrix  N  lie  in  the  interval  (-23i,2(3i)  along  the  imaginary  axis. 
If  X  is  an  eigenvalue  of  A,  it  is  contained  in  the  rectangular  region 

Re(\)  e  (0,8)  , 
Im(\)  e  (-23,23)  . 

The  Tchebychef  algorithm  was  tested  against  the  two  competing 
methods  for  values  of  3  ranging  from  .1  to  40  (see  Table  6.2). 
Figures  6.2(a)  -6.10(a)  show  the  hull  of  the  eigenvalue  estimates 
computed  by  the  Tchebychef  algorithm  and  the  best  ellipse  enclosing 
the  approximate  hull.   For  small  3  the  matrix  A  is  nearly  symmetric 
which  yields  a  hull  close  to  the  real  axis  (see  Figure  6.2(a)).   For 
large  3  the  matrix  A  is  nearly  antisymmetric  which  yields  a  nearly 
verticle  spectrum  (see  Figure  6.10(a)). 

Table  6.2 


Figure  # 

(3 

J.XJ.J- 

1 

d 

OJ.O.X 

1 
c 

■d 

-iiiaj. 

1 
C 

Convergence 
Factor 

Rate  of 
Congergence 

ICYC 

6.2 

.1 

4.0 

3.87 

1+.69 

k.66 

.9075 

.04215 

20 

6.3 

.4 

4.0 

3.87 

If.  01 

3.87 

.9502 

.02217 

40 

6.4 

.8 

k.O 

0 

3.9^ 

3.31 

'  -9737 

.01155 

20 

6.5 

2 

k.O 

0 

3.93 

3.09i 

.9571 

.01906 

20 

6.6 

4 

k.O 

0 

4.05 

8.86i 

.9558 

.01964 

20 

6-7 

8 

k.O 

15i 

3.85 

l8.44i 

•  9324 

.03039 

20 

6.8 

10 

k.O 

l4.l4i 

4.05 

20.l8i 

.9494 

.02254 

20 

6.9 

20 

k.O 

31.62i 

k.2k 

40.391 

•  9595 

.02024 

20 

6.10 

4o 

k.O 

75.00i 

5.28 

85.281 

.9769 

.01015 

20 

2 
\i 

0 
-ll 

-2._ 
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Figure  6.2(a)   Best  ellipse  for  3  =  .1 
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Figure  6.2(b)  Work  required  for  p3  =  .1 
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Figure  6.3(a)   Best  ellipse  for  p  =  .K 
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Figure  6.3(b)      Work  required  for  [3  =    .k 


Figure  6.U(a)   Best  ellipse  for  p  =  .8 
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Figure  6.U(b)  Work  required  for  p  =  .8 
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Figure  6.5(a)   Best  ellipse  for  f3  =  2 
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Figure  6.5(b)  Work  required  for  p  =  2 


Figure  6.6(a)   Best  ellipse  for  3  =  h 
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Figure  6.6(b)     Work  required   for  |3  =  k 
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Figure  6.7(a)   Best  ellipse  for  3  = 
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Figure  6.7(b)  Work  required  for  p  =  8 
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Figure  6.8(a)      Best  ellipse   for  (3  =   10 
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Figure  6.8(b)     Work  required  for  3  =  10 
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Figure  6.9(a)   Best  ellipse  for  3  =  20 
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Figure  6.9(b)     Work  required  for  p  =  20 
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Figure  6.10(a)   Best  ellipse  for  3  =  hO 
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Figure  6.10(b)  Work  required  for  ^>  =  k-0 
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Figures  6.2(b)  -6.10(b)  show  a  comparison  of  the  error 
suppression  of  the  two  methods,  Tchebychef  and  Bidiagonalization, 
versus  the  work  required.  Error  suppression  was  measured  by  the  log 


of  the  relative  error;  that  is,  by  tog.(||e  ||/||eJ|)  where  ||e  ||  is  the  I   - 

th 
norm  of  the  error  vector  at  the  n   step.  Work  was  measured  in  terms 

of  the  number  of  Tchebychef  iterations.   Recall  that  Bidiagonalization 

required  about  wice  as  much  work  per  iterative  step. 

The  initial  choice  of  parameters  d  and  c2  was  based  on  the 

rectangle  known  to  contain  the  spectrum  of  A.   In  some  cases  a  rather 

poor  choice  was  used  to  show  that  the  adaptive  procedure  would  still 

work.   Table  6.2  contains  the  initial  parameters  as  well  as  the  computed 

o 
optimal  parameters  for  each  test.   (Recall  that  c2  =  c  .)   The  convergence 

factor  associated  with  the  best  ellipse  and  rate  of  convergence  are  also 

given  in  Table  6.2.   (The  rate  of  convergence  is  -tog.  of  the  convergence 

factor. ) 

T 
For  A  nearly  summetric,  the  condition  of  A  A  is  significantly 

worse  than  the  condition  of  A.   Since  Bidiagonalization  is  sensitive  to 

m 

the  condition  of  A  A,  it  did  rather  poorly  for  nearly  symmetric  A  (see 
Figures  6.2(b)  -6.4(b)).   For  large  p  the  two  methods  were  more 
comparable,  but  the  Tchebychef  method  still  held  a  significant  advantage 
(see  Figure  6.10(b)).   For  these  tests  the  mesh  space  used  was  h  =  1. 
The  symmetric  part  of  A  is  multiplied  by  a  factor  of  l/h  ,  and  the 
antisymmetric  part,  associated  with  the  first  order  terms,  is  multiplied 
by  a  factor  of  l/h.  For  smaller  h  the  matrix  would  be  more  nearly 
symmetric,  the  type  of  system  for  which  the  Tchebychef  algorithm  holds 
the  greatest  advantage. 


In  some  tests  the  region  of  convergence  associated  with 
the  initial  choice  of  parameters  d  and  c2  did  not  include  the  entire 
spectrum.  The  error  grew  in  the  first  cycle.  Eigenvalue  estimates 
were  extracted,  however,  and  the  solution  vector  was  set  back  to  the 
initial  guess.  No  progress  was  made  toward  the  solution,  but  better 
values  for  d  and  c2  were  found  (see  Figures  6.6,  6.8,  6.9).   In  the 
test  with  3  =  k   (see  Figure  6.6)  two  cycles  were  required  before  enough 
information  was  obtained  to  choose  a  d  and  c2  that  produced  convergence. 

In  the  tests  in  which  the  region  of  convergence  associated 
with  the  initial  choice  of  parameters  d  and  c2  contained  the  entire 
spectrum,  the  very  rapid  initial  drop  in  error  was  due  to  the  suppression 
of  eigenvalues  near  the  center  of  the  region.  The  asymptotic  rate  was 
achieved  after  this  initial  drop  (see  Figures  6.2,  6.3,  6.k,   6.5,  6.7). 

b)   Singular  Problems 

Suppose  the  matrix  A  has  eigenvalues  in  the  open  right  half 
plane  except  for  an  eigenvalue  at  X  =   0.   The  Tchebychef  algorithm  can 
be  applied  to  this  case.   Since  the  scaled  and  translated  Tchebychef 
polynomials  are  all  normalized  at  the  origin  (P  (0)  =  l),  the  component 
of  the  initial  guess  in  the  direction  of  the  null  space  of  the  operator  A 
will  not  be  altered  by  the  Tchebychef  algorithm.   Indeed,  if  the  parameters 
d  and  c2  are  chosen  to  fit  the  convex  hull  of  nonzero  eigenvalues,  then  a 
solution  will  be  found  containing  the  component  of  the  null  space  present 
in  the  initial  guess.   The  error  in  the  direction  of  the  nonzero  eigen- 
vectors will  be  suppressed  as  usual. 

A  problem  would  arise,  however,  if  the  adaptive  procedure  were 
to  yield  an  approximation  to  the  zero  eigenvalue .   In  this  case  the  value 
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nearly  zero  would  be  admitted  to  the  approximate  hull  and  cause  a 
choice  of  d  and  c2  that  would  yield  a  very  slow  rate  of  convergence. 
However,  the  eigenvalue  estimates  are  based  upon  the  successive 
residuals,  and  if  the  target  vector  b  is  in  the  range  of  A,  then  the 
residuals  are  orthogonal  to  the  null  space  of  the  matrix  A.   Since  the 
residuals  contain  no  eigenvectors  associated  with  the  zero  eigenvalue, 
the  adaptive  procedure  will  not  yield  an  approximation  to  the  zero 
eigenvalue. 

To  test  the  Tchebychef  algorithm  on  a  singular  system, 
consider  the  differential  operator  above  with  Neuman  boundary  conditions. 
If  Q,(x, y)  =  0,  then  the  constant  function  is  a  solution  of  the  homogeneous 
problem 


(Mx,y)  T-)  -  t-  (Mx>y)  r-)  +  Bi(x>y)  r- +  B?(x>y)  f-W>y>  = 

1  dx         dy       d  ay  L  dx         e  dy/ 


If  we  approximate  this  operator  on  the  same  grid  as  in  part  a), 
the  matrix  A  will  be  the  same  except  that  the  diagonal  must  be  altered 
so  that  the  constant  vector  is  in  the  null  space  of  A.   As  before  let 

A1(x,y)  =  A2(x,y)  =  1, 
Bx(x,y)  =  B2(x,y)  =  (3  . 

Tests  were  run  with  3=1  and  3  =  10.   Figures  6.11(a)  and 
6.12(a)  show  the  approximate  hull  of  the  nonzero  eigenvalues  and  the 
best  ellipse  enclosing  it.   Figures  6.11(b)  and  6.12(b)  show  the  error 
suppression  of  the  two  methods,  Tchebychef  and  Bidiagonalization,  versus 
the  work  required.  Error  was  measured  orthogonal  to  the  null  space  of  A. 


Figure  6.11(a)   Best  ellipse  for  3=1 
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Figure  6.11(b)  Work  required  for  3=1 
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Figure  6.12(a)  Best  ellipse  for  p  =  10 
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Figure  6.12(b)  Work  required  for  $  =  10 


Starting  values  for  the  parameters  are  listed  in  fable  I  .'  aa  ./<rll 
as  the  computed  optimal  parameters.   The  convergence  factor  and  rate 
of  convergence  also  appear  in  Table  6.3. 


Table 

6.3 

Figure  # 

P 

Initial 
'd     c1 

Optimal 
•d       c1 

Convergence 
Factor 

Rate  of 
Convergence 

icyc: 

11 

1 

k.O 

0 

3-93 

3.12 

.98+2 

. OO69I 

20 

12 

10 

k.o 

i6i 

if. 27 

19. Hi 

•9592 

.01809 

20 

c)   Factorization  Methods 

A  recent  innovation  in  iterative  methods  is  factorization 
(Stone  [20];  Dupont,  Kendall,  and  Rachford  [5];  Dupont  [6];  Kim  [1+]; 
Bracha  [2]).   For  the  system  Ax  =  b,  the  object  is  to  find  a  matrix  B  such 
that  (A+B)   is  easily  computable  and  the  equivalent  system 
(A+B)  Ax  =  (A+B)  b  has  improved  condition.   Little  is  known  about 
the  spectrum  of  (A+B)~  A,  making  an  adaptive  procedure  attractive. 

One  factorization,  called  the  Strongly  Implicit  Procedure  (SIP), 
is  applicable  to  any  5-point  difference  matrix  (Stone  [20]).   In  this 
procedure  a  singular  matrix  B(a),  dependent  on  the  parameter  a,    is  chosen 
so  that  (A+B (a))  =  L  •  U,  where  L  is  lower  triangular  and  U  is  upper 
triangular.   The  equivalent  system  is 

(u'1L"1A)x  =  IT1!,"^  . 


13k 

When  used  in  conjunction  with  SIP,  the  Tchebychef  algorithm 

is  computationally  unchanged  except  for  computing  the  residual. 

Without  SIP  the  residual  is  given  by  r  =  b  -Ax  .  With  SIP  the  residual 

J  n       n 

becomes  r  =  U  L  (b  -Ax  ).   This  requires  ION  multiplications  and  ION 
additions  as  opposed  to  5N.  Thus,  Tchebychef  plus  SIP  requires  12. JE 
multiplications  and  13-7N  additions  per  step  as  opposed  to  7-7N  and 
8.7N. 

Notice  that 

(A+B(a))_1A  =  (A'V+Bta)))"1  =  (i  +A~1B(a))"1  , 

so  that  if  id  is   an  eigenvalue  of  A     B(a),    then 

1+jU 

is  an  eigenvalue  of  (A+B(a))~  A.   If  the  eigenvalues  of  A~  B(a)  are 
bunched  around  p   =  0  then  the  eigenvalues  of  (A+B(a))~  A  will  be  bunched 
around  X  =   1.   Since  B(a)  is  singular,  ju  =  0  is  an  eigenvalue  of  A~  B(a), 
so  X  =   1  is  an  eigenvalue  of  (A+B(a))~  A. 

SIP  was  used  in  conjunction  with  the  adaptive  Tchebychef 
algorithm  on  the  matrices  from  part  a)  with  £3  =  .k   and  p  =  k   for  several 
values  of  a.   Figures  6.13(a)  and  6.1^ (a)  show  the  best  ellipse  enclosing 
the  approximate  hull  of  (A+B(a))~  A  for  the  indicated  values  of  a- 
Figures  6.13(b)  and  6.1^(b)  show  the  error  suppression  versus  the 
number  of  iterations  for  the  indicated  values  of  a.   The  iteration 
parameters  d  and  c2  were  initialized  to  correspond  to  a  circle  centered 
at  X  =   1;  that  is,  d  =  1,  c2  =  0. 

Notice  that  the  eigenvalues  of  A  are  transformed  into  eigen- 
values of  (A+B(a))~  A  that  are  bunched  about  X  =  1.      Compare  the  ellipses 
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of  Figures  6.13(a)  and  6.14(a)  with  the  ellipses  of  Figures  6.3(a) 
and  6.6(a).  It  is  clear  from  Figures  6.13(b)  and  6.14(b)  that  the 
conditions  of  the  systems  have  been  greatly  improved.  Compare  the 
rates  of  convergence  for  the  same  matrices  in  Table  6.2  and  Table  6.4 


Figure  #     3  a  Factor        Convergence  ICYCLE 

13  .    .4  .3  -7897  .1025  20 

13  .4  .7  .6185  .2086  20 

13  .4  1.0  .4371  -359^  20 
Ik  k  .3  .2479  .6057  20 

14  4  .7  .5371  -2699  20 
14  4  l.O  .8319  .0799      "  20 


Little  is  known  about  the  effect  of  the  parameter  a  upon  the 
spectrum  of  (A+B(a))~  A  for  nonsymmetric  A  (for  symmetric  A  see  Dupont 
[5])-  It  is  apparent  from  the  results  that  the  condition  of  the  matrix 
(A+B(a))~  A  is  very  susceptible  to  the  choice  of  a.   Notice  how  the 
ellipses  in  Figures  6.13(a)  and  6.14(a)  differ  with  different  a.      The 
accuracy  of  these  ellipses  is  questionable,  however,  as  convergence 
occurred  too  quickly  to  gain  adequate  information  about  the  eigenvalue 
structure . 

The  greatly  improved  condition  and  rapid  convergence  achieved 
with  SIP  more  than  justifies  the  extra  work  per  step  and  presents  the 
explanation  of  the  role  of  a  as  an  interesting  open  question. 
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d)   Nonhomogeneous  Region 

In  order  to  test  the  algorithm  on  a  more  difficult  system,  the 
region  0  <  x  <  k±,    0  <  y  <  kl   was  broken  up  into  subregions  (see 
Figure  6.15).   The  functions  A1(x,y),  Ag(x,y),  B1(x,y),  B  (x,y)  we: 
given  values  as  follows: 
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Region  i 

\ 

A1(x,y)  = 

10 

A2(x,y)  = 

10 

B1(x,y)  = 

0 

B2(x,y)  = 

0 

Region 

B 

A]_(x,y)  = 

10 

A2(x,y)  = 

10 

B1(x,y)  = 

10 
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10 
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A]_(x,y)  = 
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0 
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A1(x,y)  = 
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A2(x,y)  = 

100 

Bx(x,y)  = 

0 

B2(x,y)  - 

10 

A  grid  with  mesh  space  h  =  1  and  Dirichlet  boundary  conditions 
were  used  to  generate  the  associated  5-point  difference  operator,  A,  of 
dimension  l600.   The  adaptive  Tchebychef  algorithm  was  used  on  this 
system,  both  by  itself  and  with  SIP.   The  results  are  shown  in 
Figure  6.l6.   Figure  6.16(b)  shows  the  approximate  hull  of  the  matrix  A. 
Figure  6.16(a)  shows  the  approximate  hull  of  the  matrix  (A+B(a))~  A  with 
a  =  -5-   Figure  6.16(c)  shows  the  error  suppression  of  the  Tchebychef 
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Figure  6.16(b)   Spectrum  of  A 
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method,  alone  and  with  SIP,  versus  work.  Hie  work  is  measured  in  the 
number  of  Tchebychef  iterations;  the  number  of  iterations  wit]  is 
multiplied  by  a  factor  of  12/7- 

Notice  that  the  entire  spectrum  of  A  is  transformed  into  a 
spectrum  closely  bunched  about  \  =  1.   The  condition  of  the  system  is 
greatly  improved.   The  convergence  factors  and  rates  of  convergence  of 
the  two  systems  are  given  in  Table  6.5« 


• 

Table  6.5 

Initial 
■d            c1 

Optimal 
*&                  c1 

Convergence 
Factor 

Rate  of 
Convergence 

ICYCLE 

TCHEB 

ko 

0 

201.61    198.11 

.998^ 

.00070 

20 

TCHEB  +SIP 

1 

0 

.7072      .6509 

.78I+O 

.IO56 

20 

In  a  problem  such  as  this,  as  much  is  known  about  the  spectrum 
of  (A+B(a))~  A  as  is  known  about  the  spectrum  of  the  matrix  A.   The 
initial  choice  of  parameters,  as  well  as  the  computed  optimal  choice 
of  parameters,  is  given  in  Table  6.5.   In  both  cases  the  adaptive 
procedure  was  able  to  recover  from  poor  initial  parameters  and  produce 
good  estimates  of  the  optimal  parameters. 

The  choice  of  a  =  -5  was  quite  arbitrary.  With  the  value 
a  =   1.0,  the  matrix  (A+B(a))~  A  had  an  eigenvalue  with  negative  real  part, 
The  improvement  in  the  condition  of  the  system  with  a  =  .5  is  significant 
enough  to  warrant  further  study  into  the  role  of  a  in  this  factorization. 


ll+2 

6.k     Summary 

The  adaptive  Tchebchef  algorithm  has  qualities  that  promote 
its  use  on  nonsymmetric  systems  whose  eigenvalues  lie  in  the  right 
half  plane.   The  method  does  not  depend  upon  any  special  structure  of 
the  matrix.   It  can  be  used  on  finite  element  matrices  and  difference 
matrices,  as  well  as  in  conjunction  with  factorization  methods. 

The  Tchebychef  algorithm  requires  less  than  half  of  the  work 

per  iterative  step  that  is  required  by  the  two  competing  methods: 

T 
Bidiagonalization  and  Conjugate  Gradients  on  A  A.   The  additional 

storage  required  by  the  adaptive  procedure  of  the  Tchebychef  algorithm 

becomes  less  of  a  factor  as  the  number  of  diagonals  needed  to  store  the 

matrix  A  increases. 

In  the  tests  that  were  run,  the  Tchebychef  algorithm  was 

considerably  faster  than  the  competing  methods.   This  is  due,  in  part, 

to  the  fact  that  it  is  sensitive  to  the  condition  of  A,  whereas  the 

T 
competing  methods  are  sensitive  to  the  condition  of  A  A.   This  advantage 

was  especially  apparent  on  nearly  symmetric  systems. 

Another  quality  of  the  Tchebychef  method  is  stability.   The 

Tchebychef  algorithm  makes  steady,  predictable  progress  toward  the 

solution.  At  each  step  the  error  is  multiplied  by  a  factor  that  is  no 

greater  than  the  convergence  factor.   This  factor  is  determined  by  the 

choice  of  parameters  d  and  c2  and  the  spectrum  of  the  matrix.  Any  round-off 

error  will  be  multiplied  by  this  factor  on  the  next  step.   The  competing 

methods,  on  the  other  hand,  depend  upon  orthogonality  and  A-orthogonality 

T 
(A  A-orthogonality) .   The  rate  of  error  reduction  of  these  methods  is 

unpredictable  at  best,  and  error  reduction  may  become  very  slow  if 

orthogonality  breaks  down  (Hestenes  and  Stiefel  [12]). 


The  adaptive  procedure  was  shown  to  be  able  to  produce  go 
estimates  of  the  optimal  iteration  parameters  with  virtually  no  a 
priori  knowledge  of  the  spectrum  of  the  matrix  A.  Although  the 
stability  of  the  adaptive  procedure  could  not  be  guaranteed  (Section  5-1), 
in  all  of  the  tests  the  procedure  behaved  as  expected.  All  of  the 
eigenvalue  estimates  lay  inside  the  true  hull  of  the  eigenvalues,  and, 
after  good  parameters  were  found,  further  eigenvalue  estimates  lay 
near  the  center  of  the  approximate  hull  as  predicted.   This  may  be  due 
in  part  to  the  absence  of  nonlinear  elementary  divisors  of  large 
dimension  in  the  test  matrices.   Further  refinement  of  the  adaptive 
procedure  may  be  necessary  in  the  presence  of  nonlinear  elementary 
divisors  of  large  dimension. 

In  addition  to  finding  good  iteration  parameters,  the  adaptive 
procedure  provides  information  about  the  spectrum  of  the  matrix  that  may 
be  useful  to  the  user.  One  such  use  is  in  the  choice  of  the  error  bound 
(Section  6.1) . 

The  adaptive  Tchebychef  algorithm  showed  great  promise  when 
used  in  conjunction  with  SIP.  Although  the  properties  of  the 
parameter  a  are  not  understood,  the  significant  improvement  in  the 
condition  of  the  systems  tested  showed  the  potential  of  this  factorization 
method . 

The  Tchebychef  algorithm  could  be  used  in  conjunction  with 
other  factorization  methods  as  well.  Y.  J.  Kim  [Ik]   has  shown  that  some 
finite  element  matrices  can  be  factored  in  much  the  same  way  as  the 
5-point  difference  operators  are  factored  by  SIP.   It  has  also  been 
suggested  that  the  Tchebychef  algorithm  could  be  used  in  conjunction 


with  Symmetric  Successive  Over  Relaxation.   In  each  of  these 
factorizations,  the  adaptive  property  makes  the  Tchebychef  method 
attractive  because  little  is  known  about  the  spectrum  of  the  equivalent 
system. 
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C  SUBRGUTTNE  TCHE B < A, X , B, R , DX, XLA ST f P LA  ST t S tCH , 

C  1                             ERBND,  ICYCLE,  IMAX) 

C  THIS  SUBROUTINE  SOLVES  THE  SYSTEM  AX=B,  RETURNING  THE      * 

C  SOtUTICN,  X.  ^HE  PCSITIVE  HULL  CF  tj-e  EIGENVALUES  OF  A  IS  * 

C  DETERMINED  DYNAMICALLY  AND  RETURNED  IN  CH.  * 

c * 

C  * 

C  THE  USE*  MUST  SUPPLY  THE  SLBPCLTINE  NSYMAX (Y  ,X  ,  A  )  ,           * 

C  WHICH  PREFORMS  THE  MATRIX  VECTGD  MULT  I PLL IC ATIGN :           * 

C  Y  -  AX.                                          * 

C  * 

C  THE  INPUT  PARAMETERS  APE:                                       * 

C  * 

C  A(N,NC)           THE  MATRIX                                  * 

C  X(N)              THE  INITIAL  GEUSS                         * 

C  B(N)              THE  TARGET  VECTOR                         * 

C  * 

C  EPENC                              THE    ACCEPTABLE    EFROR                                                       * 

C  ICYCLE                            THE    NUNBER    CF    ITERATIVE    STEPS    EETWEEN         * 

C  ADAPTIVE    PPCCEDURES                                                         * 

C  IT^AX                                THE    MAXIMUM    NUMBER    OF     IThRATINS    ALLOWED    * 

C  * 

C  CCMMON  /APEAl/  N,NC  :                                        * 

C  N        THE  DIMENSION                               * 

C  NC       THE  NUMBER  CF  COLUMNS  NEECEC  TC  STCRE    * 

C  THE  MATRIX  A(N,NC)                         * 

C  * 

C  CCMMON  /APEA2/  D,C2,A2,FC  :                                 * 

C  D        THE  CENTER  CF  THE  ELLIPSE                 * 

C  C2       THE  FCCAl  LENGTH  SCUAPED                  * 

C  * 

c * 

C  * 

C  OTFER  VARIABLES  A«E:                                         * 

C  * 

C  P  (N)              THE  PESIDUAL                      * 

C  DXCNJ             THE  CHANGE  IN  X  AT  EACH  STEP     * 

C  XLAST(N)          X  AT  THE  START  OF  THE  CYCLE      * 

C  RLAST(N)          P  AT  THE  START  OF  THE  CYCLE      * 

C  S<N,4)            THE  LAST  A  RESIDUALS  CF          * 

C  EACH  CYCLE                         * 

C  * 

C  CH(25,2)          EIGENVALUE  ESTIMATES             * 

C  TCH(2i>,2)        LINK  LISTING  OF  CH                * 

C  * 

C  A2                PEAL  AXIS  OF  THE  BEST  ELLIPSE    * 

C  CF                CCNVERGENCE  FACTOR  OF  THE        * 


c 
c 
c 
c 
c 
C 
c 
c 


c 
c 
c 
c 

c 
c 

c 
c 


c 
c 
c 
c 


RQ 

P  1 

n2 


PEST     ELL  IPSf 

FEMDU:    NCPM    C*     RfcSICUAL 

RESIOU    AT     THk     STAPT    OF 

THE    CYCLt 

IltPATlUN     FAPAKETtPS 


* 
* 
* 
ft 

* 

* 

SURCCUT INE    TCFFB(A,X,R,P,DX,XLAST,rLAST,S,CH, 
1  EFBNC,  ICYCLE,  I"AX ) 

TMPIICIT    REAL*fi( A-HfC-Z ) 

fC^CN    /A9EAL/    N,ND 

COKNCN     /AREA?/    D,C2,A2,FC 

COMMON    /AP6A3/     I  CM , I F I o ST , I F o EE 

CINEf^  TCN    A(NlKOIlX(N)fE(M|ft(N)lDX(Nh 
1  XL AST(N)  ,PLAST(N)  ,S(N,4) 

DI^ENSICN    C"(25,2) ,  ICM 25, 2) 

BEGIN 

INITIAL  IZE    Cl-t  ICI-.lFlPSTf  IFPEE,A2tFC 

CALL    TNTT(CH) 

INITIALIZE     ISTEP 
!STE°    =    0 

CALCULATE    ~FE     INITIAL    RESIDUAL 
6  CALL    N^YMAX(C ,X ,A) 

DO    8    1=1, N 

MI)     =    B(I)     -    R(I) 
C  CCNTTNUF 

OUTPUT    THE     INITIAL       PESIDU 

CALL     r.L"fPLT(P,PSD,IS">  LP) 

SAVE     TNI~IAL    X    ANC    R 
PC    =    PSO 
CALL    LASTX( X,P , PSO, > LAST, PL  AST, PC) 

BEGIN    STTEFfcl     ITEFATICN 

TNITTALTZF    OX    AND    P  AP  AMtT  Cr-  S    PI    AND    P2 
LC  DC     15     I     =     1  ,N 

DX(I )     =    "(I ) /D 
lb  CCNTTNUE 


Dl    =     2.D0/D 
F2     =    O.CO 
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c 

C  *AI N    LCCP 

20  DO    40    I    =     l.ICYCLE 

C 

C  STCRE    LAST    FCUP    PESIOUALS     IN    S 

IFd.GT.ICYCLE-4)    GC    TO    25 
GO    TO    3C 
C  THEN    STCRE    R 

25  IC0L  =  ICYCLE-H-1 

CO  26  K=l,N 

S(Kf ICCL)=P(K» 

26  CONTINUE 
C 

C  CC^PUTE  NEW  X 

30  00  32  J=1,N 

XCJ)  =  XU  )  ♦  0X(J  ) 
32  CONTINUE 

C 
C  CGPFUTE  NEW  P 

CALL  NSYMAX(P,X,A) 
DO  34  J=1,N 

R(J)  =  B(J  )  -  F(J  ) 
34  CONTINUE 

C 
C  COMPUTE  NEW  PARANETEPS 

P2  =  C2*P1/(4.D0*D  -  C2*Pl) 
PI  ^=  (l.CO  ♦  P2)/D 
C 
C  COMPUTE  NEW  OX 

OC  36  J=1,N 

CX(  J)  ■  P1*R(J  )  ♦  P2*CX(J) 
36  CONTINUE 

C 
C  UPDATE  ISTEP 


C 

C  OUTPUT  RESIDU 


ISTEP  =  ISTEP+1 

:  RESIDU 

CALL    OUTPUT(P,R SO, ISTEP) 


40  CCNTJNUE 

C 

C  END    NAIN    LCCF 

C 

C  END    STIEFEL 

C 

C  TEST    ^C   HALT 

TF(PSD.LT.EPBNQ)    GO    TO    <30 
IF(  IrTEF.GE.IMAX)    GC    TO    <3C 


c 
c 

f. 
c 
c 


c 

c 
c 


<*b 


5C 


pl-GIN    ADAPTIVE     FRi.CECLPE 

CALL     ADA°t<s  ,fl  ,CH,  I  FLAG  I 

TF    D    C°     C2    HAVE     PEFN    CHANGLC    RLSTAFT    STIbFEU 
CThtPWTSE    RESUME    STIbFFL 

TF  (  IFLAG.FCC.l  )    GO    TO    4b 

GO    T[j    5C 
THEN    <<ECTART    STiEFEL 

CALL    LASTX(X,R,RSC,XIAST,PLAST,FC) 

GO    T  Q    I  C 
EL  CF    RESUME    STIEFEL 
GC    TO    2  0 

END    ACAP7  IVF    PPOCECtJR  E 


9  C    P  E  T  L'ft  N 
ENC 


SLBRLLTTNE  CLTPLT (R , P S C , IP'  E P ) 

C    ^HtS  SUBROUTINE  GU7PUTS  !>E  ITERATION  PARAMETERS  C  ANO  C2,   * 
C    THE  EXPECTED  CONVERGENCE  FACTCP,  FC,  ANC  TFE  NCPM  OF  P.      * 

Q   ********#*********************************»>:*****************#** 

IHPLTCT"  REAL*8( A-HtC-2l 
CC^CN  /APEAl/  N,ND 
CCHNCN  /ARFA2/  C,C2,A2tFC 
DINFNSJCN  2<N) 


BE^TN 

CALL     INPTMP  ,f  ,B«C,N) 
P«D    =    DSC DT (PSD) 


IOC 


ENC 


WPITEC6,  LOO)     IS"rEP,C,C2,FC,PSD 
FOP^AT(»     ISTEP     =,,I3,3X,'D    =    • , D12 . b, 3X , • C2 
D12.5,3X,'FC    =     •  ,D12. 5t3X,'RESICL 


-  I 


=    •  t0l2.b) 


RETURN 
END 
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SUBRCIFIN 

(^  ************** 


c 

ThIS    SUBROU 

c 

CF    X    AND    P 

c 

************* 

IMPLICIT 

CCFFCN    /A 

DIMENSION 

c 

c 

BEGU 

IF(PS 

THEN 


10 


15 


20 


ELSE 


25 

C  END 

C 

30    RETU&N 
END 


E    LASTX(X,R,RSC, XL  AST, RL AST, RO  ) 

**********************************  ****** 

TINE       PEPLACES    X    AND    R    BY    PREVIOLS    VALLE 
IF    THE    PESICU,    PSD,     HAS    GPCWN, 

**************************************** 

PEAL*8( A-H,0-Z) 
REA1/    N,ND 
X(N),P(N) ,XLAST(N) ,RLAST(N) 


D.GT.PC)    GC    TC    10 

GO    TQ    20 

CHANGE    X    ANC    P 

DO    15    I=1,N 

XU  )=XLAST(  I) 
R(I >=PLAST(T) 

CONTTNUE 

PSC=PO 

GC   tc    3  0 

UPDATE     XLAST    AND    PLAST 

PO    =    RSC 

DC    25    1=1, N 

XLASTU  )=X(U 
RL AST(I )=P( I) 

CONTINUE 


***** 

S        * 

* 
***** 


c 
c 
c 
c 


c 

C 


SUBPCUTINE     INFPC(X,Y,  FPCCN) 

********************************************************** 
THTS    S'JBPOUTINE    FINDS    THE    INNEF    PRODUCT    CF    THE    TWO         * 

VECTCP^:     PRCD    =    <X(N),Y(N)>  * 

********************************************************** 

IMPLICIT  REAL*8( A-H,C-Z) 
DIMENSICN  X(N) ,Y(N) 


BEGIN 

FRCD  =  O.CO 
DO  10  1=1, N 
PPCD  = 
1C      CCNTINLE 
RETUPN 
ENC 


PPOO  +  X( I)*Y( I ) 


SlRPCUriNE     TMT<o) 

i  ThIS    SUBROUTINE    INITIALIZE*    CH , t ch, I F I P ST, I F PE E  ,  A2  ,»  C    * 

C    THi    HULL     IS    CIV  FN     r>~E    VALUES     D*Ci     C-C    ANC    TEE     LLLIPM      IS     * 
C    THE    LINE     BETWEEN    THEM.  * 

INFLICT"'    PEAL*8 ( A-H,C-Z) 
r.OKMCN     /AREA2/    C,C2,A2,FC 
Cf/VPCN    /APEA3/     ICF,  IFIPST,  IFREE 
OIMENSICN    CH<25,2) ,ICH(25,2) 

BECIN 


C 

c 

c 
c 


c 
c 


INITIALIZE    TCH 

DO    5    1=1,25 

IfH(I,  I) 

5  TCH(T  ,2) 


I-L 
H-i 


10 


20 


INITIALIZE    CF,TFIPST, IFPEE, A2tFC 
IFIRST    =    1 
!F(C2.LE.0.0DC)    GC    TQ     10 

^C    TC    20 
THEN 

IF°EE    =    2 

ICM1  ,1)    =    1 

!CH(l,2)     =    1 

CHI  1,1)    =    D 

CHIi.2)    =    CSQRK-C2) 

A2    =    C.COC 

CC    TO    30 
ELSE    (C2.GT.0) 

IFPEE     =    3 


TCF(ltl)    = 
TCH  (I,?)     = 
ICH(2, 1)    - 
ICH(2,2)    = 
CH(1,1)     = 
CH(  1,2)    = 
CF(2,1 )     = 
CH(2,2)     = 
A2    =    C2 


1 

2 

1 

2 
0-DSCPTCC2) 
0.000 

D40SCRT (C2) 
C.ODC 


30    FC     =    C.COC 

PETUPN 
END 
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SL'BRCL 

C**#******  ** 

C  THIS    S 

C  BASED  ON  T 
C  VALUES,  CE 
C    MATPEX    S. 

£*********** 

IMPLTC 
CCMKIN 

COMMON 
DTMENS 
OIMEN^ 


C 
C 

C 
C 

C 
C 


c 

c 
c 

c 
c 
c 

c 
c 


TINE  ADAPT(S,R ,CH, 

****************** 

UBPOUTINE  FINDS  Th 
HE  CCNVEX  HULL  CF 
NEGATED  FROM  THE  R 

****************** 
TT  REAL*8( A-H,0-Z) 

/AREA1/  N,ND 

/APEA2/  C,C2,A2,F 
ICN  S(N,4) ,R(N) 
TCN  CH(25y2)tQ<4), 


IFLAG) 

******************************** 

E    BEST    ITERATION    PARAMETERS  * 

A    SET    CF    APPRCXIMATE    EIGEN-  * 

ESIDUAL    VECTORS    STORED    IN    THE      * 

******************************** 


10 


FINC    THE    LEAST    SQUARE    SO 
CALL    LSTSC(S,C,R) 


EV<4,2) 

LUTION    TO       SC*R=C 


FINC   TFE    EIGENVALUE    APPROXIMATIONS 
CALL    EVS(C,EV,IRT) 

IF    NO    NEW    EIGENVALUES    RESUME    STIEFEL    ITERATION 
IF(  IPT.EC.O)    GO   TC    5 

GC    TO    IG 
THEN 

IFLAG    =    0 
GC    TO    2C 


ACC    THE    NEW    APPPCXIMATIC 
AND    FORM    THE    NEW    CONVEX 
CALL    HULL <CH,EV,IFLA 

IF     THE    EIGENVALUE    AFPPCX 

THFN    RESUME    STIEFEL     ITER 

IF(IFLAG.EC.O)  GC  TC 


NS  TO  THE  PPEVIUUS  CONVEX  HULL 

HULL 

G,IFT) 

IMATICNS  ACC  NO  NEW  INFCPMATICN 
ATION 
20 


FINC  the 
CALL 


BEST  PARAMETERS 
ELLI P<CH) 


FOR  THE  HULL 


20 


RETURN 
END 


SIHPCLMNE    LSTSCCStQfM 

C  THIS    SUBPOJ"MhF    FINDS    Tf-t    LEAST    SQUARES    SOLUTION    OF       * 

C  THE  SYSTFM,  $*Q  =  -B,  BY  SCLVINC  THE  SYSTEM,  STS*C  =  B,  * 
C  WHERE  ST  =  S-TPANSPUSE  AN  C  e  =-ST*R.  A  BI 01  AGU,M  AL I  I  AT  IUN  * 
C    ROUTINF     I  *    LSFD     TC    SCLVE    THE    SMALL    SYSTEM.  * 

£********444********************  *******■■#*  a****************** 

TMFLTCIT    PEAL*8 ( A-H,C-Z ) 

CCMHCN     /APFAl/    N,ND 

OIMFNMON    S(N,4),Q(4),P(N),STS(4,4),B(4) 


C 
C 


C 

c 


c 
c 
c 
c 
c 
c 


10 


CCMPLTE    Sf*S 
CO    1C    1=1,4 

CC    10    .1  =  1,4 

CALL    INPFC(S  (1,1  )  ,SU,  J),STS(I,J),N) 
STS( J, I  )    =    STS(t  ,J) 
CCN^  TMJE 


COMPLTE    R 

CC    20    1=1,4 

CALL    TNPCC(SUtI  )  ,  ^  ,  B  (  I)  ,N) 
8(T  )    =-e(  I  ) 

20  CCNTINUE 

NCPMALTZE    THE     SYSTEM 
ALPHA     =    S^Sd  ,1) 
OC    3C    !=1,4 

6(1)    =     e(I)/ALPHA 
DC    30    -1  =  1,4 

STS<I,J)     =    STSU  ,J)/ALFHA 
30  CCNTINUE 

SOLVE    THF     SYSTEM    STS*C    =    B 

***ANY    LTBPARY    POUTINE    MAY    BF    USED    TO  SOLVE    *** 

***               the    4X4    *YSTEM     :       STS*C    =    B.  *** 

***                BICTAG    APPEARS     JN    APPENDIX    B.  *** 

CALL    RIDTAG (STS ,C,E) 


FE^UPN 
END 
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SUBROUTINE    EVS(G,EV,IPT) 

C  THIS    SUBPCUTINE    FINDS    THE    RCCTS    OF    THE    POLYNOMIAL  * 

C    WITH    THE    COEFICIENTS:     1,QC  11 »Q( 2) , QI31 tC<4! .    THESE    PCCTS         * 
C    APE    THEN    TRANSFORMED    INTO    EIGENVALUES    OF    THE    MATREX    A.  * 

TMPLICTT    REAL*8< A-H,Q-Z) 

CCNKCN  /ARE*2/  C,C2,A2,FC 

DIMENSTCN    AA( 51,3) ,RCCTS(2,5U) ,CD<51> 

DIPENSTCN    EV ( A, 2) , Q ( 4  ) 

BEGIN 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 

c 
c 


FINC  PCCTS 

***ANY  LIBRARY  ROUTINE  MAY  BE  USED  TO  FIND  THE  PCCTS. 
***7HE  PCCTS  SHOULD  BE  STUPED  IN  EV<4,2)  FROM  THE  FRONT. 
***FV(I,1)  IS  THE  FEAL  PART  AND  EV(I,2)  THE  INAGINAFY  PART 
***TFE  VARIABLE  I»T  INDICATES  THE  NUMER  CF  &CCTS  FCUND. 
***  PCSR  APPEARS  IN  APPENDIX  C. 


10 


AA(  1,  1)  =  l.CDO 
DC  10  T=l ,4 

AA(H-1,1)  = 

CONTINUE 
IOEG  =  4 


C(  I  ) 


14 


16 


CALL    RSSR(AA,P00TS»IDEG,1C, 15, 1. CD-4, 1. CD-6 ,DDI 

IPT    =    4-IDEG 
CO    14    T=1,IPT 

FV(  1,1  )    =    RCCTS(l,5-i ) 

EV( I ,2)     =    POCTS(2,5-I) 
CCNTINUE 


TRANSFORM    THE    ROOTS    TO    EIGENVALUES    CF    A 
G    =    D+DSCPT(C*D-C2) 
K=IRT 
H  =  0 
CC    20    T=i,K 

T1  =  EV(I  ,l)*EV(I,i)+EV(I,2)*EV(I,2» 
IF(Tl.LT.CABS(C2)/(G*G) )    GO    TC    16 
GC    tc    13 

THEN  DISCARD  RCCT 
IRT=IFT-! 

GC  TC  20 
ELCE  FIND  EIGENVALUE  ESTIMATE 


c 
c 


le 


20 


22 


24 


30 
ICC 
200 
30C 

40 


M  =T  I  ♦  l 

FVUI,l)  =  C-.5CO*EVU,l)*(C2/(Tl*<")*G) 
EV(II  ,2I»-.5CC*EV( I,2)*(C2/(Ti*G)-G) 
CONTINUE 

OITPLT    NFU    FIGENVALUE     ESTIMATE^ 
TF<  IFT.FC.O)    CO    TO    22 

GC    TC    24 
THEN 

WTTE(6,  300) 

GC    TC   40 
El<E 

wnJTE  (o,  100) 
nc    30    I  =1  ,I&T 

WRITE(6,200)     EV(I  ,1),EV(I  ,2) 
CCNTIMJE 
F0"WA^'0THE    NEW    EIGENVALUE    ESTIMATES     APE:') 
FC°MAT(«     «,5X,  012.5,*     ♦     I*     SO  12. 5) 
FC»*A7(»0NG    NEW     EIGENVALUE     ESTIMATES') 
PETLTN 
ENC 
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C*** 

c 

C  TH 
C  IF 
C  HU 


C 

c 
c 
c 
c 

c 

c 


SUBROUTINE  HULL(CF,EV,  IFLAG,  IPT  ) 

THIS  SUBROUTINE  FORMS  THE  POSITIVE  HULL  CF  THE  UNION   CF  * 

E    PREVIOUS    POSITIVE    HULL    ANC   THE    NEW    EIGENVALUE    ESTIMATES.    * 

NONE  OF  THE  NEW  EIGENVALLE  ESTIMATES  ARE  IN  THE  POSITIVE   * 

LL,     IFLAG    IS    RETURNED    AS    0.  * 

IMPLICIT    REAL*8( A-H.O-Z) 
CO*CN    /ARE/2/    C,C2,A2,FC 
COMMCN    /AREA3/    I CH, I FI P ST , I  FREE 
DIMENSION    CM25,2),ICM25,2),EV(4,2) 

BEGIN 


SET     IFLAG 

IFLAG    =    0 

PLACF    NEW    EIGENVALUE    ESTIMATES     IN   THE    HJl L 
TF(  IRT.E0.01    GO    TO    <5C 
CC    50    I=i,IRT 

IF(FV(I,2).LT.0.0D0)     GC   TO    50 


C 

C 

C 
C 


TEST=(D-EV(  It  II  )*<D-EV<  1 ,  1 )  )  ♦  <  A2-C2  )♦ 

EV(I ,2)*EV(I ,2)*A2-A2*(A2-C2) 
IF( TEST.LT.O.COCJ    GC    TC    2 

GC   TC    4 
THEN    EV(I,-)     IS    NCT     IN    THEN    HULL 

GO    TO    50 
ELSF    SET     IFLAG 

IFLAG    =    I 

PUT    EV    IN    NEXT    CPEN    SPOT 
J=IFREE 

CH< J,1)=EV(I ,1) 
CM  J,2)=EV(I,2) 
IFPEE=    ICH(J,2) 

LINK     NEW    MEMBEP     IN    CRCER 

IFCCH(Jfl).LE.CMIFIPST,m    GO    TO    8 

GO   TO    1C 
rj-EN 

ICH(Jtl)     =    J 
ICH(J,2)    =    IFJRST 
ICHUFIRST.l)    =    J 
IFIPST    =    J 
GO    TO    5C 


10 
15 


20 


K     =     T  F I P°l 

K     =     ICH(  K  ,2) 

IF(CH(J,  1)  ,LT.Ch(  Kf  1)  )    GC    TC    2C 

GC  TL  2  5 
THFN    LINK 

ICH(J,  1  I     =     ICF(K,  1) 

IO(Jt2)     =     K 

ICH(  TCH(K,i  )  ,2)     =     J 

ICKK,  1 |     =     J 

GC    TC    50 


C 

C 

C 
C 

c 
c 


25 


30 


50 


e  ND    OF 


CCNTINUE 


TF( I CF  (Kf 2)  .FC.K  ) 

GC    TC    15 
THEN    LINK    CN 

ICH(K ,2) 
ICH( J,i) 

ICHC J»2) 

GC   TC    50 
LINK 


GC    TU     30 


END 

=  J 
=  K 
=     J 


IF     NONE    OF    THE    NEfc    EV'S    WEF  E     PLACED     IN    TFE    HULL    CET|J3N 
TF( TFLAG.EQ.O)    GO    TC    90 


FCFH    NFW    HULL 

K    =     ICF(TFIR5T,2) 


6C 


65 


7C 


IF( TCH( K,2) .E.C.K)     GC    TC    9C 

TEST= (CHK,2)-0(  ICMK,  1)  ,  2)  )*<CHUCMKf  2)  t  l)-CM(K,l)  ) 
-(Ch(IC(-(K,2),2  )-CH(K,2  )  )  *<  CF  (K  ,  1  )-CH  (  I  CF  (K,  1  ) ,  1 )  > 


TF( 


GC    TO    65 


!CH(K ,2) 
ICF(K, 1) 


TEST. LE.C. CDC) 

GC  ~0  7C 
THEN    LNLINK     CH(Kf-J 

ICH(TCH( K,IJ,2)    = 

TCH(TCH(K,2),  1)     = 

TCH( Kf2)     =    IFREE 

TF^EE    =    K 

K    =     !CH(K,1 ) 

!F(ICH(K,  D.EC.K)     K    =    ICH(K,2) 

HC  '0  60 
EL^F     *CVE    tc     NEXT    K 

K=     JCHlKf2) 

GO    to    60 


9C    PE^L'SN 
END 


l6o 


SUBPCUTINE       ELLIP(CH) 

C  THIS    SUBROUTINE    FINOS    THE    OPTIMAL    ITERATION    PARAMETERS    * 

C    FOP    THE    FCSITTVE    HULL    CH.    IT    RETAINS    ONLY    The    KEY    ELEMENTS    * 
C     IN    THE    HLLL    AND    OLTPLTS    THESE,  * 

C++****** 4 **+ 4+ 4*** *+*4+**X*+1 +*++*+++++*+++*+++++++++++++++++ 

IMPLICIT    REAL*8< A-H,C-Z) 
COMMON    /AREA2/    0,C2fA2,FC 
CCM^CN    /AREA3/     ICF, I F  IRST , IFREE 
DI  ME  f^SI  CN    CH  (25,2)  ,ICH(  25,2  )  ,ICCL(  25) 
C 

C  EEGIN 

C 

C  FIRST    SEE    TF    ANY    PAIPfclSE    BEST    IS    THE    CPTIHUP 

J    =     TFI&ST 
10  IF( ICH( J, 2) .EC. J  )    GC    TO    30 

K    =    J 

15  XFUCMKt21.EQ.KJ    GO    TO    25 

K    =    ICH(K,2) 

CALL  TWOPT(CH,J,K,IFLAG) 
IF(IFLAG.EO.l)  GO  TO  15 
C  TFFN  TWOPT  FAILEC 

C 

C  TEST  TO  SEE  IF  THIS  PAIPWISE 

C  BEST  IS  CPTIML. 

I  =  IFIPST 
C 

It  IFCI.EC.J  .CP.  I.EC.K)  GL  TO  18 

GC  TC  2C 
C  tj-EN  SKIP  TEST 

16  GC  TC  22 
C 

2  0  TEST=(C-CH(  1,1)  )*<D-CH<  I  ,  1)  )  * ( A2-C2)  ♦ 

1  CMI*2)*CMI92)*A2-42*(42-C2) 

IFUEST.GT.l.CD-8)  GC  TC  21 
GC  TC  22 
C  THEN  TEST  FAILS 

21  GO  TC  15 
C 

22  IF(  ICH( I  ,2). EC. I)  GC  TC  23 

GO  TC  24 
C  THEN  TEST   WCFKS 

23  WPITE(6,5CC) 

GC  TO  90 

C 

C  UPDATE  I 

24  I  =  ICH(I,2) 

GC  TC  16 


L6l 


c 
c 


c 
c 

c 
c 
c 
c 


c 
c 

c 
c 


END    LF     TEST 
2  5  J    =     I  C  M  J  ,  2  ) 

gc   Tn    10 

END    PATPWISE    SCARCE 

FINC    TME     BEST    TH&EE    VALLE     PCI  N^ 

INITIALIZE     TCCL     AND    TFC 
30  DO     22    !  =1,25 

TCUL ( I)    =    0 
3?  CCNTIME 

TFC    =     l.OC 

TPY  FACE  'F^EE  VALUE  FCTNT 
j  =  IF'^1 

35      TFUCE(  TCH(J,2),2)  .EQ.!CH(J  ,  2)  )  GO  TU  7C 

K  =  TCH(J,2) 
40  TF  UCH<  K,2).EC.K)  GC  TC  65 

I.  =  K 

45  TF(ICH<L,2).EC.l)  GC  TC  60 

L  =  ICH(L,2) 

CALL  TEFEPT (Ct,J,K,L,  IFLAG ) 
IF( IFLAG. EC. 1 )  GC  TC  45 
THEN  THeEPT  FAILED 

TEST  TC  SEE  IF  THIS  TEFEE  VALUE  PCINT 
IS  A  CANDIDATE 
I  =  IFISST 


46 

48 
50 


1 


51 
52 

53 


IFU.EQ.J  .CP.  I.EO.K  .CF.  I.E3.L)  GCTC  48 

GC  ~C  50 

THEN  SKIP  TEST 

GO  TO  52 

TES1  =  (C-CF(  1,1  )  )*<C-CH(  1,1 )  )*(A2-C2) 
♦CH(  I  ,2)*CH( I ,2)*A2-A2*(A2-C2> 

IF(TEST.GT.l.0C-8 )  GC  TC  51 
GC  TC  52 

TEEN  TEST  FAILS 
GC  TC  45 

IF(  ICH(  I, 2)  .EO. I  )  GO  TO  53 

GC  TC  58 
THEN  TEST  WCRKS 

ICOL(J)  =  I 

ICCL(K)  =  1 
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ICCL(L)    =    i 
IF(FC.LT.TFC)    GCTO    54 
GOTO    56 
C  THEN    THIS  POINT  BEST  SO  FAR 

54  TFC    =    FC 

TC    =    D 
TC2    =    C2 
TA2    =    A2 
56  GC    TO    45 

G 
C  UPDATE    I 

58  I    =     ICF(I,2) 

GC    TC    46 
C  END    OF    TEST 

C 

6C  K    =    TCH(K,2) 

GC  TO  40 
65      J  =  ICH(J,2) 
GO  TO  35 
C      ENC  OF  SEADCF 
C 

C      THE  BEST  THREE  VALUE  FCINT  FIT  HAS  BEEN  FCUNC 
70      C  =  TC 

C?  =  TC2 
A2  =  TA2 
FC  =  TFC 
C 

C      SAVE  ONLY  KEY  ELEMENTS 
I  =  IFIPST 
8C      K  =  7CHU  ,2  ) 

IF(  ICGLd  I. ECO)  GO  TO  82 
GO  TO  84 
C  THEN  LNLINK  CH(I,-) 

82  !CH( TCH(  1,11  ,2)  "    ICH(I,2) 

ICHI  ICHI I,2)tl)  =  ICFU, 1) 
ICH(I  ,2)  =  IFPEE 
TF°EE  =  I 
84      IFCK.EC.T )  GC  TC  SO 
T  =  K 
GC  TO  80 
C 

C      OL^PLT  NEW  PARAMETERS  ANC  POSITIVE  HULL 
90      WOTTE(6,600) 

WFTTE(6,7C0)  D,C2»A2,FC 
VPT^Et  6,8C0) 
I  =  IFIRST 
95  WPTTE(6,900)  CH ( I , I ) , CH ( I ,2 ) 

IF(  ICH(I  ,2). EC. I )  GC  TO  99 


c 

c 


T     =     ICH(I 
GC    to    93 


2) 


500 
60C 
700 

80C 
900 


FORMAT^ : 

FCKA'C 
FCPVAT( 

FC^MAK 

012. 5, 3X, 

Ft)»MAT(  t     THF 


THE 

THE 
C    = 


BtST  .  •  ) 


3X,'C2    =     ■  ,012.5, 3>, 'Ai    =     • 


CPT  IML     IS     A    PAIFWlSfc 

NEW  PAPA^FTEFS  ARE : • ) 

•,D12.5, 

'FC    =     ', C12.5) 

POSITIVE    HILL:* ) 


FC»MTC     •,  3X,  012.5,  '     ♦I*',D12.5) 


9S 


P  F  "  I  P  N 
ENC 
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SUBROUTINE    TWOPT(CH, J,K,TFLAG) 

C  THIS    SLBPCUTINE    FINDS    THE    OPTIMUM    ITERATION    PARAMETERS    * 

C    FOP    THE    TWO    EIGENVALUES:    CH(J,-)     ANC    CH(K,-)  * 

IMPLICIT    REAL*8< A-H,0-Z> 
COMMON    /APEA2/    C,C2tA2,FC 
CTKEf^SICN    CM25,2) 


C 
C 


C 
C 
C 


c 
c 


c 
c 


SET    THE  CONSTANTS 

/*    =  (CH(Ktl)-CH(J,  1  ))/2.0CC 

B    =  (CH<K,1)+CH(  J,l.)  )/2.0C0 

S    =  (CH(K,2I-CH(J,2))/2.CCC 

T    =  (CH(K,2)+CH(J,2))/2.0C0 

ELIMINATE    DEGENERATE    CASES 
CASE    I:    T    C*ALL 

TF(T/B.l.T.i.OD-3)    GC    TO    5 
GC    TO    10 
.    THFN 
5  D    =    B 

C2    =    A*A 

A2    =    C2 

FC    =    A/(D«-0SCPT(D*D-C2)  ) 

GC    TO    50 

CASE    II:    A    SMALL 
10  JF< A/T.LT.l.OD-3)    GC    TO    15 

GC    70    20 
THEN 
15  D    =     B 

Y    =    T+DABS(S) 

C2    =    -Y*Y 

A2    =    A*< A+DSCPT(A*A+4.C0*Y*Y) )/2.0DC 

FC    -    <DSCPT< A2)*CSCFT (A2-C2)  )/<C4CSQRT(C*C-C2)  ) 

GC    TO    50 

CASE    III!     S    SPALL 
20  TF(CABS(S)/A.LT.1.0C-3)    GC    TO    25 

GC    ^c    30 
THEN 
25  CALL     ACUBIC(A,B,T) 

GC    TC    50 

CASE     TV:    GENERAL    CASE 
30  CALL    F!FTF(A,B,S,T, IFLAC) 


50    RETURN 

END 


c 
c 

c 
c 


SlRHriTINr.  TMPFPT(CH,J,Kfl,lFl/\C) 

THIS    SUBROUTINE    FINDS    ThE     ITERATION    PARAMETERS  + 

ASSOCIATED     HTM    THE     LMCUt     EILIFSE    THRCUGH    THE    ThPH  * 

EIGENVALUES:     C1-"  U  ,- I  tCH(  K  ,- )  ,     AND    CH(L,-)  * 

IMPLICIT    9L.  AL*6(  A-H,  0-Z) 
CCMMCN    /APEA2/     C,C2,A2.FC 
DIMENSION    CH(2b,2) 

SET     TFLAC 

TFLAG    =    0 


COMPUTE    D 

U     =     CHCJf2  )*CM  J, 2  )*(CMK,  l)*CH(K,  IJ-CM  L,  1MCHIL,  I)  I 

1  tCHI L,2)*CH( Lt2)*(CH( J, 1)*CH( J,1)-CH(K,L )*CH(K, L ) ) 

2  ♦CMK,2)*CMKf  2)*(CMLt  l)*CMLf  1)-CH(J,1)*CH(  J,  II  ) 


Z    =    CH( J,2)*CH( J,2)*(CH( K,1)-CH(L ,1) ) 

1  ♦  CH(Lt2)*CML,2)*(CHIJ»  1»-CH(K,  1)) 

2  ♦CH(K,2)*CH(K,2)*(Ch(Li 1)-CH( J, I ) ) 


C 

C 


IF  (  Z  .L E.O.OCO  )    GJ    TO     10 

nc  Tc  2c 

THEN    NU    ELLIPSE    FITS    THESE    TH^EE    VALLES 
10  IF  LAC-    =     1 

GC    TC    90 

20  C    -    W/(2  .0C0*Z) 

TF( O.LE.C.CDC )    GO    TC    30 

GU    TO    40 
THFN    FLLI°<F    CONTAINS    THE    ORIGIN 
3C  I  FLAG    =     1 

GC     TO    9  0 

CCMPL'TE    OTHER     ITERATION    PARAMETERS 
4  0  A2=(CH{ J,2)*CH(J,2)*CF(L, l)*CH(K, I )*(CH(L, 1)-CH(K, 1) ) 

1  +CH( L»2i*CH<L,2)*CM J, l)*CF(K,i)*(CH(K,i)-CH(J, 1)) 

2  ♦CH(Kf 2)*CMK, 2)*CF( J, 1)*CF(L, 1)*<CH{ J, l)-CH(Ltl))  ) 

3  fl    ♦    C*D 


50 


TF(A2.GE.0*DI    GO    TG    5C 

GC    tc    60 
THEN    ELIPSE    CONTAINS    THE    CPIGIN 

IFLAH    =     l 

GC    TC    90 


60  C2=A2*<1.0C-Z/(CH< J , 1 ) *CH < L , 1 ) * < CF ( J  ,  1 > -CH ( L  ,  1  )  ) 

1  *CF(L,  1)*CH(K,  1I*(CE(L ,1I-CH(K,  II) 

2  +CH<  J, I )*CH<K, l)*(CF(K,ll-Ch(J,i)|)J 

FC=(CSC°T(/>2I+CSQCTU2-C2i  ) / ( D+OCQPT <D*D-C2 ) ) 


90    ^etlpn 
ENC 
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SU8PGLT 

Q******** **** 

C  THIS    SU 

C    WHEN    BCTH    E 
C    BEST    VALUE 
C    NCMIAL. 
£************ 

IMPIICT 
CC^NCN 


C 
C 


INE  ACUBIC(A,B,T ) 

***************** 

BPOUTINE  FINOS  TH 
IGENVALHES  HAVE  T 
CF  A2  IS  FCUND  A5 

***************** 
T  REAL*8( A-H,C-Z) 
/AREA2/  C,C2tA2.F 


******************************* 

E  BEST  ITERATION  PARAMETERS    * 

FE  SAKE  IMAGINARY  PART.  THE    * 

THE  POCT  CF  A  CUBIC  FCLY-     * 

* 
******************************* 


TO  CCMPUTE  D,A2,C2,FC 
X  =  (B 


<B*B*T*T  )*(B*B*T*T)+A*A* 
(B*B+T*T)*( B*B+T*T)*(B*B+T*T)) 
(e*E+T*T  l*(B*B+T*T J ) 


*B*A*A*A*A*T*T)*< 
(T*T-8*B)  )/( 2. DO* 
Z=(B*B*A*A*A*A*T*T \f  ( 

>1=X*DSCRT(X*X+Z*Z*Z) 
Y2=X-0S0P"'(X*X*Z*Z*Z  ) 
Y=DSTGMOABS  (Yl  )**(l  .C0/3.C0  ),Yl  )  +  CS  IGN(  CABS  <  Y2  )  ** 

(l.DO/1. DC) ,Y2) 
A2=Y-M8*e*A*A  )/( B*B*T 
C2=(A2*( A2-T*T-A*A|)/ 
FC    =    (OSOPT(  A2H-DSCR 


*T) 

(A2-A*A) 

T(A2-C2))/(C+CSCFT(C*C-C2)) 


OETLRN 

END 


SUBCClTJNfc     FTf-TMAf  H,S,Tf  IF  LAG  J 

THI<:    SUePOJTTNE    FINDS    THE    OPTIMAL     !TL»A1ION    PA°AMFTtRS  * 

FOP    TWC    CC^PLEX     EIGENVALUES.    The     EEST    VALLF    UF     L     IS     FOUND  * 

AS    THE     PCOT    OF     A    FIFTH    OEGREE     FCLYNCMlAL     MTF    COEFICrENTS  * 

PliP<!fP3iP^P5tP6.     (COMMON     /AREA*/).  * 

IMPLICIT    *FAL*8< A-H,C-Z) 

CCMMCN    /ARE*2'     E,C2,A2»FC 

CC^CN     /APEA4/    PI ,P2 tP3 tP4,P5  ,P6 


CGMPUT6    ".>E     FTFTF    DECREE    PCLYNCMIAL 
"1     =    -A*7/<: 
02    =    -A*S/T 

cj     =       $*~/A 
*4    =    -9 


Pl*fil*<»l«,2.D0*F2*P3-4.D0*R4J+P2*IP2+R3-4.0C*P41 

1  +  P4*(4.D0*r.4-2.C0*P3) 

P2  *P  1  *<  4. DO  *£4*<P3-«?4)«-P2*(12.D0*P  4- 5. DC*P  3-4. 00**2) 

1  ♦?l*(2.C0*F4-33-4.C0*P2) ) 

2  +  c2*(4.D0*P4*(R3-F4)-»-P21M2.D0*F4-P3)) 

3  -r3*r^*-r4 

P3=Pl*(R2*( C4*  <2.CC*P4-4.CG*P3 )+P2*( 3.D0*P  3-4.D0*R  4)  ) 

1  +  P  l*(P4*34*-P2*(4.CC*P2«-3.D0*&3-4.D0*P4)-Pl*R3)  ) 

2  -f  P2+P2*  (R^*c  4-P  2*P3  ) 

C4sRl*P2*o**(Rl*<3.CC*Pi-4.C0*P4)*P2*<3«CG*'»2-4.D0*R4) 
1  +  2.Q0*R4*P4) 

P5=3.C0*C l*P 1*R2*R2*F2*(2.CC*F4-P  I-R2  ) 
P  £  =  Q  i*F  l*P2*P2*P3*(Rl*P2-F4*f-4) 

CCMPLTE    Y     AND    Z 

V    =     <T-S)*(T-S)*  (E* A)*(B  +  A)-(T4S  )*  <T*S ) *< P.-A)*(B-A  ) 
Z    =     (T-S  )*(T-i;)*{B  +  A  J    -     (T*S)*("r+S)*(6-A) 


IF( S.LT.C.OCO)    GC    TO    10 

C-C    T0    20 
THEN 
IC  El    =    DMAXK-A,  (Y/(2.CC*Z)-E)  ) 

E2    =    0.00 
GC    ^C    50 
ELSE 
20  IF(Y.C"r.2.C0*E*Z  )     CO    TO    30 

GC    "*0    4  0 
TFEN 
30  El    =    O.ODO 

E2    =    A 
GC    TH     50 
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C         ELSE 

40        El  =  O.CO 

E2  =  0*INMAtY/(2.C0*Z)-e) 

GO  TO  50 
C 

C      FIND  ROOT 

C  ***ANY  LIBRARY  ROUTINE  MAV  BE  USED  TC  FIND  THE  REAL  ROOT*** 
C  ***  CF  the  FTFTH  DEGREE  PCLY  IN  THE  INTERVAL  (E1,E2).  *** 
C   ***  THE  ROOT  SHCDLC  BE  RETURNED  AS  El.  *** 

C   ***  ZEROIN  APPEARS  IN  APPENDIX  C.  *** 

C 

50     CALL  ZERC!N(E1,E2,IFLAG) 
C 

IHIFLAG.GE.3)  GC  TC  60 
GO  TC  70 
C         then  ZEROIN  FAILED 
60         IFLAG  =  1 
GC  TC  80 
C  ELCE 

7C        IFLAG  =  0 
C 

C  CCMPUTE    ITERATION    PARAMETERS 

D    =    El     ♦    B 
A2=    (E 1-P1)*(F1-P2) 
C2=    A2*(E1-P3)/(EI  ) 

FC=    (DS0RT(A2  ) +DSCFT ( A2-C2  ) ) / < D  +  DSCRT <C*C-C2  )  ) 
C 

80    FE'UFN 
END 


APPENDIX   B 
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c 
c 

c 

c 


c 
c 
c 
c 

c 
c 
c 


c 
c 

c 
c 


SUBPCLTINE    BTOI AG(A,X,B) 

THIS    SUBROUTINE    FINDS    THE    LEAST    SQUARE    SOLUTION    OF    THE       * 
SYSTEM;    A(N,M)*Xm     -    ti(N)     -    0.  * 

IMPLTCI"r    PEAl*8(*-H,C-Z) 

DIMEhSICN  A<4,4),  X ( 4 ) , E ( 4)  , U ( 4  ,6 )  ,  V ( 4  ,6 ) f VW (4  )  ,T i (4 ) ,T 2 ( 4) 

CI^EhSICN  P(6) 

N  =  A 

M  =  4 

IMAX  =  4 

BEGIN 

INITIALIZE    ISTEP 
TSTEP    =    0 

INITIALIZE    VECTQPS    ANC    CONSTANTS 
SET    x(M»  ,  VM  K)     =    0 
00    2    1  =  1, M 

X(  I  )    =    O.DO 
VU(I)     =    O.DO 
2  CONTINUE 

SET    Bl*L(N)     =    B(M 

CALL    INPftOIBtBtBlvM 

Bl    =    OSQRT(Bl) 
00    6    I =  1,N 

U(I  ,1)     =    fi(T )/Bl 
6  CCNTINUE 

INITIALIZE    A1,W,Z,T 
Al     =    O.DO 
W    =    O.DO 
Z    =    -l.CO 
T    =    l.DO/Bl 

SF^     EPNC 

EBND    =    Bl*l.GD-lO 

FIFST  ST£P:  COMPUTE  FIRST  V  AND  Al 
CALL  TMLL(V(l  ,1)  ,U( i,l)  ,A,N,M) 
CALL  INPROCVIltlltVCltlltAltMl 
*1  =  CSCRT(Al) 
DC  8  1  =  1, N 

V(T,  1)  =  V( I, 1I/A1 
8        CCNTTWJE 


IV 


c 

c 

c 

c 


c 
c 


c 

c 


^ATN    10CP 

UPDATE     I^TFP 
IC  I S7FP    ■    ISTFF*1 

CC^FUTF     NEW    X 

L     =    -Z*ei/Al 
DO     12    I  =  1,M 

X( I )    »    X(I )     f    Z*V(  I, ISTEP  ) 
12  CONTINUE 

CCNPLTE     NEW    V* 

'«     =     <T-R1*W)/A1 
CO     14     T=l,M 

VMI1     =    VMII     ♦    W*V<  I,  ISTEP) 
14  C1NTTMLE 

CC^FLTE    NEW    U    *NC     El 

CALL    A^LL(T1  ,V(l ,1  STEP) ,A,N,^) 
CG     L6     1  =  1,  N 

U(IfISTEP*ll    =    TIUJ    -    A1*U  (  I,  ISTEP  ) 
16  CONTINUE 

CALL     !NPFU(U(i,ISTEP  +  l),U<l,ISTEP4-l),Bl,N) 

si    =  nscFTien 

TEST    Bl 

IF<81.L7.EBNC)     GC    ^C    18 

GO    TG    2C 
ThEN    SCLUTICN     IS     EXACTLY    X 
16  GC    TG    ICC 

ELSE    PE-GPTFOGONALUt 
20  K    =    ISTEP+l 

DC    22    T  =  1,P 

U(  I,K  )    =    U(  I,K)/81 
2?  CCNTINUF 

OC    24    1  =  1, ISTEP 

CALL     IhPPOCU  UVK)  VU(  It  It  VP<  I  )tM) 
24  CONTINUE 

OC    26    1  =  1, ISTEP 
CC    26    J=1,N 

U(  J,K)=U(  J,  K)-FU  )*U(J.  I  ) 
26  CONTINUE 

CALL     INPFC (U Cl.KJ ,U( 1,K  ),B2,M) 
82    =    0SCPT1B2) 
DC    28     '  =  1,M 

U  (I,K1     =    U(  I,K  )/E2 
28  CONTINUE 

UPDA'E    7 
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T    =    -7*AI/B1 
C 
C  TEST    ISTEP 

TFCTSTEP.GE.IPJX)    C-Q   TO    100 
C 
C  COMPUTE    NEW    V    AND    Al 

CALL    TPUL(T2,UUVISTEP  +  1),A,NVM) 
DO    3C    1*1, M 

V(I,IS7EP+i)    =    T2U)    -    B1*V(I  , ISTEP) 
3C  CCNTINLE 

CALL    INPFOCVU  ,ISTEPU)  ,Vt  1,  ISTEP*  I  >  ,  Alt  PI 
Al    =    DSCFTUl) 
C 
C  TEST    Al 

IFUULT.EBNC    .OF,    ISTEP.  EC.  4)    GO    TO    32 
GO    TO    38 
C  THEN    SOLUTICN    IS    X=X-G*VW 

32  G    =    B1*Z/(B1*W-T) 

OC    34    1*1 f> 

X(  I  )    =    X(I  )    -    G*VMI1 
34  CCNTINUE 

GO    TU    ICC 
C  EL  CE    RE-ORTHOGONALIZE 

38  K    =    ISTFP+1 

DC    40    I  «i,M 

V(  I,K  )    =    V( J,K)/A1 
40  CCNTINUE 

DO    42    1*1, ISTEP 

CALL     TNPR0(V(1,K),V(  1,I),P(I),M) 
42  CCNTINUE 

DC    44    1*1, ISTEP 
DC    44    J=1,M 

V(J,K)=V<  J,K)-F(I  )*V<J,  I) 
44  CONTINUE 

CALL    INFFC(V(1,K),V(  l,K  ),A2,M) 
A2    =    DSCPTU2) 
CO    46    1*1, M 

V  (I,K)    =    V(I,K)/A2 
46  CONTINUE 

C 
C  FEPEAT    LCCP 

GO    TO    1C 
C  ENC    CF    LOOP 

C 

100    RETLPN 
ENC 


L73 


c 
c 
c 
c 


<tRoriTINE     AMU  JY,X,A,N,M 

******* | 4**4*4********* **** *** ******** ******************** 

THI<  SUBPCUTINF  LLHS  TFE  MATRIX  MULTIPLICATION:        • 
Y(N)   =  A( N,M)*X( M)  .  * 

*****4-*4*******************4****4********, **************** 

TKPLICTT  =*EAL*8  (  *-H,  f-Z  ) 
CIMENSICN  A( N,M)  ,X( M)  ,Y(N) 


10 

20 


BEGIN 
OC 


CCN 

PETLPN 
ENC 


2J  I  =1  ,N 

Y( I )  =  C.DO 
DC  1 C  K  =  l  ,M 
Y(  I)  = 
CONTINUE 

TTNUE 


VCIJ  ♦  ACT  ,K)*X(  K) 


c 

c 

c 
c 
c 


c 
c 


SUBTCUTINF  ^^UL (Y,X ,A, N,M 

***********  *********************************************** 

THIS    *'JBPOUTINE    COES    THF    MATREX    MULTIPLICATION:  * 

Ym    =     AT(MfM*X(N),  * 

WHERE    AT     IS    THF    TPANSFCSE    OF    A.  * 

******  4  **4***4**********>,'*»M  *************  **********  ******* 

TMFLTCI7    0FAL*8 ( A-H,C-Z ) 
DIMENSION    AC N »M) ,X( N 1 1 VI  HI 


10 
20 


BEGIN 
DC 


CCN 

RETUPN 
ENC 


20    K  =  l,M 

Y(K)    =     0.00 
UL    10    1=1»N 

Y(K)     =     Y(K)     ♦    A(I  ,K)*X(  I) 

CONTINUE 
T  I  MJ  F 


I7h 


APPENDIX  C 


MlflFOUTINF     FSSP     (A,  "  GOT  S  ,  DEGF  EF  ,  M  tP  b  AX  ,  CI  l_T  A ,  F  PS  I  L        ,C) 

DIMENSION     A(5LI,         T  A(  5  L  ,  3  )  ,  P  OHT  S  (  ^  ,  50  J  ,  0  (  5  I)  ,P  C.MUD  (  5  C  )  , 
IN0NCT(25I ,MNONP~ (25) 

PEAI   *8  A ,POOT<tD,FOMOD,DELTA,FFSlL 

INTFGEF     CFGPEE 

NCUP -DEGREE 
'    NL=NCUR 

CALL     POG-TSQ    I  Av  lAtNCUPfM) 

CALL     FEAL^T  (  Af  I  A  ,  ¥  ,  NCUP  ,DFLT  A  ,  E°S  I  L        ,RCMCD,IA, 

1N0NPT, MNON°T, NCO,RCCTS) 

TF     (NCG.EO.O)     GO    T0     12 

CALL    CCMPRT  (A,  IA,PC*CDtPCCTS  t*SMNCNFTfNCNF""f 


1  NCU, DELTA, FPML 

IF     (NCUP .EQ.O  )    CO    *0     12 
TF     (NCLP.EO.NL)     M  =  *  +  i 
TF    (M.LE.MMAX)    GO    TC    7 
GO    T0     13 

12  CALL    FECCN    (ROOTS, A, 

13  DE^EF=NCUP 
PETUDN 

END 


,NCUP) 


CCECFEE) 


SUBFOU'-INE    ROCTSC     I  At  I  At  NCUP  ♦  MM) 

DIMENSION    A( 51,3)  ,  IA(  51, 3) 

FFAL*P  A,X 

Nl  =  NCUR  +  l 

00     1    J=i,Nl 

A( J,3)=A(J,1  ) 

IA( J, 3  »  =  C 

CALL     SCAL     (A(  J,  3),  IA(  J, 3  )  ) 

CONTINUE 

DO    9    M=ltMM 

DO    2    J=1,N1 

A( J, 2) =A( J, 3) 

IA( J,2)=TA( J, 3) 

A( J,3)=0.0 

CONTtNue 

DO    6    J=1,N1 

KM=    MTNO     (N1-J,J-1) 

IF    (KM.EO.C)    GO    TO    5 

DO    4    L=L,KM 

JL=J-L 

Jl.P  =  J  +  L 

X=A( JL,2)*A(JLP,2  ) 

LR=    MOD    (L,2) 

IF    (i.P.FQ.i)     X=-X 

IX  =  IA( JL,2  )+IA(JLP,2) 

CALL    ADD     (A( J, 3) ,IA( J, 3) ,X,IX) 

CONTINUE 
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A( J,3)=2.0*A(J,3) 

X=A( J,2)**2 

IX=IA(  J, 2)+I A(J,2) 

CALL    ADC    (A( J,3  )t IA(J,3),X, IX) 

JR  =    MOD       <J,2) 

IF     (JR.EC.O)     A<Jf 3)=-A< J,3) 

CONTINUE 

C0NTINLE 

RETURN 

END 


SUBFOUTINE    FEALRT  < A , I  A, t, NCU ", DELT A, EPS  IL       ,POMOD » MRGMGD, 

1NQNCT, MNONRT,NCUt ROOTS) 

DIMENSION    A(51t3)fIA(51t3)t ROOTS (2 , 50 ) , ROMGD( 50 ) , MFC* CD( 5  0 ) , 
IN0NFTI25) ,MN0NRT(25) ,  IPIV(5l) 

PEAL*8  A, POCTStPOMOD.T.XNtW,  Q  ,  DELTA  ,EPSI L 

NCUP1*NCUR-H 

DO    1    I*1,NCUR1 
1    IPIV(I)=1 

IF    INCUR. EQ. I )    GO    TO    8 

DO    7    I=2,NCUP 

IF     (A(I,  3)  .EO.0.0  )    GO    70    6 

T=  AUt2)**2/MIt3) 

TF     (MODd  ,2).EQ.O)    T=-T 

IT=TA(  X,2)+IA(I, 2)-IA(I,  3) 

XN=T*64.0**IT 

IF( DABSC l.-XN).LT.DEl TA)     GO    TO    7 

6  IPIVd  )=0 

7  CONTINUE 

8  1=1 
T'*=0 

DO     12    I2=2fNCURl 

IF     <IPTVU2).E0.0  )    GO    TO     12 

T4=!4+l 

I  1  =  1  2-1 

R0MCD(I4)  =  DABS( A(  I2,3)/A(  !t3)  ) 

I*30MOO  =  IA(I2,3)-IA(I,3  ) 

IF    (ROMODI I4).EQ.C.CJ    GO    TC    11 

T=2.0**H*FLCAT(I1  ) 

XN=(DLCG(R0M00U4)  )  +  DFLOAT  ( I  PCKCC  )  *4 

POMODC  I4)  =  DEXP(  XN  ) 

11  MRCyCD(I4)=Il 
T=T2 

12  CONTINUE 
Q=0.0 
NCO  =  0 

DO  22  1=  It  14 
KL=I4*l-T 

W=-POMOD(KL) 
I5=MROMOC(KL ) 
DO  20  J=1,I5 


I  58  88308335967 186C0)/T 


i  77 


i«  cali  test  (a,   w,o,ncup ,ornrn( KL)  ,FPSIL 

TF  (K  .  EQ.O)  Gf  Tf  2  1 
ROO^M I ,NCUR) =-W 
PUOTS( 2tNCUc ) =0.0 
DO  16  L=2,NCUR1 
16  A( L, 1) =A<L, i) -ACL-i  .1 >*W 

20  NCUP=NCUR-l 
GO  ^C  21 

21  W=-W 

IF     (W.GT.O.O)     GO    TC     14 
NC0  =  NCOl 
NONFT(NCO) =KL 
MN0NPT(NC01=I5*l-J 

22  CONTINUE 
PETUPN 
END 


»K) 
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SUeFOUTINE    COMPPT  ( A, I  A, ROMOD, ROOTS, M, MNCNRT,NONRT , 

i  NCO, DELTA, EPSIL       ,NC'JP) 

DIMENSION    A(5  1),  IA( 51  )  ,POMOO( 50  ) , ROOTS ( 2 , 5C)  ,SR < 48, 3) , IS* (4  8 ,3 ) , 
1SP0M0D<47) ,SRC0TS<2,47),MNCNPT(25),N0NPT<25),NS0NFT(23),MSNGPT(23! 

PEAL*8  A,P0MCD,R0CT3,SP,SPCMCD,SRCCTS,  W, C  , DELTA  ,  EPS IL 

DO    23    1=1, NCO 

JA=NCNRT(I) 

U=MN0NPT(T)/2 

DO    11    J=1,I1 

CALL    SLBRES     I AfIA *KCUR ,SR , ISRf PCMOD < JA ) I 

NSCUR=NCUR-3+4/NCUR 

J2  =  NSCUR 

IF     (NSCUR.NE.  I)    GC    ^C    9 

SROOTSUt  1)  =  -SR(2,1 )/SP<  i, i) 

NSCUR=0 

GO    T0     il 
9    CALl     PCOTSO    (SR,ISP,N^CUR,MJ 

CALL    PEALRT  ( SP  ,  ISP , ^, NS CUP ,DELT A, EPS IL       ,SPCMOD,     ISR, 

lNSONRT,MSNORT,NSCC,SPrCTS) 

11  LL=J2-NSCUR 
L=l 

12  L2=J2*1-L 

IF    (LL.EC.O)     SRLCTS(1,L2)=0.0 
!F(DABS( SROOTS< L,L2)) .GE.2.0)    GC    TU    16 
W=SF00TS(1,L2 )*RQMCC(JA) 
Q    =FCMCD( JA)*FGPLC( JA) 

CALL    TEST    <A,  W,0    ,NCUR, PCMOD < JA)  ,  EPSI L       ,K) 

TF     <K  .EQ.O)    GO    TO     16 
P00"S(1,NCUP  )  =-V</2.0 
R00^S<2,NCUR )=DSCR"M0-(W/2.0)**2) 
POOTS( 1, NCUR-l)*RCCT« <1,NCUR) 
POOTS( 2,NCUR-l) =-PCCTS (2 ,NCUP) 
CALL    OUADIV    (NCUP,A,         W,0) 
NC'JP=NCUR-2 
GO    T0    22 
16    l=L+l 

IF     (L.LE.LL)    GC    TC    12 

22  CONTINUE 

23  CONTINUE 
RETURN 
END 


SUBFOUTTNE    TEST     (A, 
DIMENSION    A(5l), 
&EAL*R  A,E,W, 

R( ?)=0.0 
P(3)=A(  I) 

DO   2     I=i ,N 

Bt  l)=B(2) 

B(2)=B<3) 

BC3»=A(I*l)-W*E(2)-C*R(il 

CONTINUE 

T|  I  >  =  0.0 

T<2)=0.0 

Nl=N*l 

E(l  )  =  PCMCD+2.0*EP$TL*PCMOD 

E(2)=P0PCD«-EP$IL*FC*Cn 

DO    7    I'1,N1 


UtQt  NfPDKCCbCSIL        ,K) 

Bl  3)  ,  7(2  )  ,F.U) 

E,        DIE, 0, R0MCU,T,f PSIL 


18 


T(l)»T(i)*E(l  HOAe$(MI)  ) 
T(  2)«T(2)*E(  2)+DA8$(  A(I  )  ) 

CONTINUE 

DIF=T(l)-T<2) 

K  =  0 

IF     (O.NE.0.0)    GO    T0    18 

TF(DIF.GE.DABS(8<3))>     K=l 

PETUPN 

!F    (DIF**2.GE  .  (  0*8  (  2  )**2«-W*B  (  2  J  *B(  3  ) +B  (  2)**2)  ) 

PF""tjRN 

END 


K=l 


i8o 


SUBROUTINE    SUBRES    ( A  ,1  A, N, SR ,1 SR , RCMCD) 

CIMENSTON    A<51),       IM51),       SRU8),        ISR(A8),       C( 51 I ,B( 48, 3) 

REAL*8  A,SR,POCD,       T,C,       E 

Nl»N+l 

T=1.0 

DO     1    I=l,Nl 

J=N1-I*1 

CU)=A(J)       *T 

T=T*ROMOD 
1    CONTTNUE 

IF     (N.LE.4)     GC    TC    16 

N2=N-2 

CO    3    1=1, N2 

BC I ,11=0.0 
3    CONTINUE 

B(1,2)=C(1 ) 

DO    8    1=2, N2 

B(  l,3)  =  C(I  >-B(l,l  ) 

DO    5    J=2,I 

B( J, 3)=-B( J-l ,2)-6( J, I) 
5    CONTINUE 

DO    7    J=1,I 

IF     (J.NE.I)    B( J,L)=B( J, 2) 

B( J,2)=B(J,3  ) 

7  CONTINUF 

8  CONTINUE 

$s(N2>       =C(N)-B(1,3) 

S«MN2-ll        =-C<M)-B<2,3) 

DO    10    J=3,N2 

K=N2+1-J 

SR<K)       =-BCJ,3) 
10    CONTINUE 

RETUPN 
16    SF (  1)       =C(  I) 

ss<2>       =-C(2) 

IF     (N.EQ.2)     RETURN 

S*(3)       =C(3)-C< 1)-C<5) 

PETURN 

END 


5UR0UITINE    RECCN     <FCL7S,A,  D,M 

DIMFNSTON    ROOTS! 2? SOI vO(  511 

PEAL*8  PCCTS,A,C,       T,U 

DO     1    1*1, N 
C(I  )*0.0 

1    CONTINUE 

D(N  +  1 )=i.O 

NS=N*l 

DU     10    I =1 ,N 

NL=NS-I 

IF     <RCCTS<2,I))    3,7,10 
3    x=ROOTS( 1,1 )**2+RCCTS<2,I )  **2 

U=2. 0*ROOTS(  1, I ) 

DO    5    J  =  M,N 

o(  j-i»=t*o<  j-n-u*n(  j)*o(  j*  i) 

5    C0NTINUE 

0( N)=T*D(N)-U*D  (N+l) 
GO    to    S 

7  T=-POOTS( Lt I ) 
00    8    J=NL,N 

D(  J  )*T*D( J)+D< J+l) 

8  CONTINUE 

S    0(N+i) =T*0(N+1) 

10  CONTTNUE 

00    li     IT=1,NS 
0(TT  I  =D( IT )*A 

11  C0NT!NUE 
PETUPN 
END 


MIBPOUTINF    OUADIV     (N,A,  W,C) 

DIMFNSION  A< 51) 

rEAL*8  A,       W,C,XN 

Nl=N+l 

DO    3    1*2, Nl 

XN=A(I-1 )*H 

IF    (I.NE.2)    XN*>N>AU-2)*Q 

A(T  )=A(I  )-XN 

CONTINUE 

^ETU^N 

ENC 


182 


11 

12 


SUBROUTINE    ACC    <X,IX,Y,IY) 

PEAL*8  X,Y,  B 

IF    (X.NE.CC)    GO    TC    3 

X  =  Y 

TX  =  TY 

GO    tq    13 

IF    (Y.EQ.O.O)    RETURN 

B  =  Y 

IDIFF=IX-IY 

TF    UDTFF)    5,12,8 

B=X 

X=Y 

IX=TY 

!OTFF=-IOIFF 

IF    (  IDIFF.GE.14)    GC    TC    13 

CO    11     I=1,IDIFF 

B=B/64.0 

CONTINUE 

X  =  X+B 


13    CALL    SCAL       (X,IX) 

PETURN 
END 


SUBROUTINE    SCAL       (X, IX ) 

PEAL* 8  X,V 

Y=DABS<X) 

IF    (  Y.f^E.O.O)    GO    TC    3 

I  X=0 

PETURN 

IF     (Y.LE.64.0)    GC    TC    5 

Y=Y/64.0 

IX=TX*1 

GO    to    3 

IF    (Y.GE. 0.015625)    GC    TO    7 

Y=Y*64.0 

JX=*X-1 
GO    to    5 

IF     (X.LT.O.O)    Y*-Y 

X=Y 

PETURN 
ENC 


L83 


5UBP0UTINE    ZERCIMe,C,  TFLAG) 

!MPITCIT    RFAL*8(A-H,0-Z) 

COMMON    /APEA4/     P  1  ,  F  2 1  P3»  P4f  °  5t  P6 

F(Z)=     (<(<Pl*Z«-P2)*Z  +  P3)*Z*P4)*Z  +  P3MZ*P6 

ABStPP     =     O.DO 

CELFRR     =    0.00 

1'=9.0D-15 

PF  =  OMAXl(t>EL  E^R  tU  ) 

IC=0 

ACB5=0ABS(B-C) 

A=C 

FA=F(A> 

FB=F(B) 

FC=FA 

K0UNT=2 

FX=HMAX1 (CABS (Fe) ,rAB^ (FC) ) 

TF(OABM FC) .GF.DABS(FB) )  GU    '0    2 

A=B 

FA=FB 

B  =  C 

FR-FC 

C  =  A 

FC=FA 

CM3=0.500*<OB) 

ACM6=DABS (CMB ) 

*0L=RE*0ABS(B)+ABcFP0 

IF{ ACMB.tE.TOL )CQ    T0     3 

IF(KCUM.GE.500)GC    TO    12 

P=( B-A)*FB 

0=FA-FB 

TF(P.GE.0.00 )G0    TG    3 

P^-P 

A=B 

FA  =  FB 

IC^IO  1 

IF(  !C.L.T.4)G0    TO    A 

TF(8.0C0*ACMB.GE.ACB<:  )CC    TC    6 

TC=0 

ACP5"  =  ACMB 

IF(o.G^.OABS(C)*TCL)GC    ^C    5 

B=B+DSIGN(TOl,CMB) 

GO    TO    7 

IF( P.GE-CMB*01GC    TC    6 

p=B+P/0 

GO   ^c    7 


i&k 


6  B  =  G.5DC*(OB) 

7  FB=F(B) 
IF(FB.EC.O.DO  )GC    TC    9 
KOlNT=KCUNT+l 

IF(DSIGN( l.ODO, FB).NF.D$IGN< l.OOO, FC))GC    TO    1 

C=A 

FC=FA 

GO    TO    1 

8  IF(DSTGN(i.ODO,FB).EQ.CSIGN( i.0CO,FC))GO  TO  li 
IF{DABS<FB) .GT.FX )GC  TC  10 

IFLAG*1 
RETURN 

9  IFLAG=2 
RETURN 

10         TFLAG=3 

RETURN 
li  IFLAG=A 

RETURN 
12  TFLAG-5 

RETURN 

END 
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